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ABSTRACT 

By assuming the validity of the Rankine -Hugoniot conservation relations for 
interplanetary type shocks in an isotropic medium it is demonstrated that im- 
proved shock normals can be calculated by employing a least squares technique 
to combined magnetic field and plasma data from a single spacecraft. The scheme 
uses only those conservation relations (six in number) which are devoid of pres- 
sure and temperature terms. Transforming these equations cast for a shock 
frame of reference into an arbitrary frame reduces the system to three indepen- 
dent "overdetermination" equations. These three equations constitute a three 
parameter redundancy among the eleven measured parameters of the system: 

Bj , B 2 , W (= V 2 - Vj ) , p x , and p 2 , where subscripts 1 and 2 refer to before and 
after the shock respectively. By exploiting this redundancy in the cases of 
simulated shocks, whose basic noiseless characteristics are known exactly, it 
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has been shown for many realistic examples, through the minimization of a least- 
squares loss function, that the normals are calculated with error improvements 
of factors of about 3 or so over calculations using the magnetic field alone. 

<C|A corrected normal and improved shock parameters are then obtained for 
a real case: the August 29, 1966 (Pioneer 7) shock. An appendix provides a 
listing of the complete computer programs used in obtaining the best estimate 
shock parameters, the shock surface normals, and the associated error cones. 

The scheme should prove useful in examining the shape of a shock surface 
whenever data for a shock event are available from two or more spacecraft widely 
separated in solar longitude. 


iv 



IMPROVED SHOCK NORMALS 


OBTAINED FROM COMBINED MAGNETIC FIELD 
AND PLASMA DATA FROM A SINGLE SPACECRAFT 

I. INTRODUCTION 

In space research there is increasing need for obtaining more accurate 
shock surface normals. This report presents a method of improving the calcula- 
tion of oblique shock normals, over previous single spacecraft methods, by using 
combined magnetic field and plasma data from a lone spacecraft. One important 
reason for obtaining improved single spacecraft shock normals, which is of 
particular interest to the authors, lies in the observational study of the shape of 
interplanetary shock surfaces. For this type of study it is presently rare to 
obtain reliable data from two spacecraft widely separated in solar longitude, a 
situation necessary for this surface shape determination, much less from three 
or more spacecraft. If N(N >3) number of spacecraft useful for this sort of 
study do exist, that is, do reliably see the same shock surface, then one can be 
reasonably sure that N - 1 of them, or at best N - 2, will be located in the near 
earth region. And in no case in the forseeable future will a situation exist where- 
by two spacecraft will be located far from the earth in solar longitude and at 
the same time remain in close proximity to each other. By close proximity we 
mean at least close enough to each other to see the same shock normal almost 
at the same time (i.e. with a time difference on the order tens or hundreds of 
minutes). Hence, we must be satisfied with reliably calculating the shock normal 
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from single spacecraft data, especially for the far-from-earth spacecraft. 

Figure 1 describes this situation where the far spacecraft shown is a Pioneer 
spacecraft (or could be considered any other solar orbiting probe), and the 
plane of the figure is approximately the ecliptic plane. The near-earth space- 
craft could represent one or more of the Explorers or any other capable space- 
craft in that region. If 6 , the difference in solar longitude of the two spacecraft, 
is sufficiently large, then for around the time of the shock sighting the two cal- 
culated shock normals n p and n £ should suffer a difference great enough to 
yield a respectable determination of the shock surface's curvature. In order 
to accomplish this the error angles associated with the estimates of the normals, 
represented in the figure as error cones of cone angles a p and a E respectively, 
should each be significantly smaller than 9 / 2 . 

Previously when one wished to calculate the shock normal from the data of 
a single spacecraft the magnetic field alone was used in the expression 

_ _ (b l x bJ x (I, - bJ 

" ' "IK x B a ) x (B 2 — B,)| ’ <I ‘ 1) 

where B x and B 2 are the magnetic fields before and after the shock respectively. 
The plus or minus sign ambiguity is clarified once the sign of the plasma density 
change is ascertained, but quantitative knowledge of the density is not required. 
Expression (I - 1) rests on the so-called coplanarity theorem (Colburn and 
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Figure 1. Two spacecraft, widely separated in solar longitude Q t each detecting different portions 
of the same shock surface- The plane of the Figure is approximately the ecliptic plane 
which contains the unit vectors R and T, orthogonal to each other and to N which is 
normal to the ecliptic- The quantities n and a refer to the unit norma I vector to the shock 
and its associated error cone angle, respectively. 



Sonett, 1966) which in principle holds exactly. However, the values for the 
vectors B t and B 2 undergo fluctuations, and therefore straightforward average 
values are often used. If these averages are markedly different from the 
"actual" magnetic field values required by Expression (I - 1), then the effect 
of the errors in the B’ swill usually be magnified by the non-linear nature of 
the expression yielding a rather unreliable estimate of the normal. This is a 
particularly serious problem when the angle between B : and B 2 is small, say 10° 
or less, and the rms deviation of the field quantities is substantial, character- 
istically say 0.6 y for the components of Bj and perhaps 1.0 y or larger for the 
components of B 2 , around the shock transition region. For this case, where 
and B 2 may be ~6 y and^lly respectively, the error cone angles for Bj and B 2 
themselves are each about as large as the average angle between them. These 
errors are then propagated by way of the two factors Bj xB 2 and B 2 -Bj which 
join to increase the error in the final calculation (I - 1), leading often to a very 
inaccurate result. Conversely, a small increase in the accuracies of B : and B 2 
should result in an even greater improvement in the accuracy of the shock 
normal’s estimate. By utilizing the associated plasma data along with physical 
relationships connecting the plasma quantities to the magnetic field quantities 
we expect to obtain at least some improvement in the estimates of Bj and B 2 . 

And this improvement, however small, will propagate its way through Equation 
(I - 1) to provide hopefully a significant improvement of the shock normal’s 
estimate. It is expected, in most eases, that this improvement will occur even 
if the plasma data is acquired with poorer accuracy than the field data. 
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In Part HE we describe this method of best fitting both plasma and magnetic 
field data which will be done by exactly satisfying the basic conservation relation- 
ships of the shock system. The best-fit magnetic field parameters are then used 
with Expression (I - 1) to obtain the so-called best estimate shock normal. Con- 
sistent with the proof of (I - 1) via the coplanarity theorem, which follows directly 
from the conservation equations, Part m uses these equations through a "best 
estimate" scheme to obtain the "proper" values to he used in (I - 1). It must he 
stressed that without such a scheme it is not at all clear what values for Bj and 
E 2 are to be used in (I -1). Surely shock parameters do not appear as simple 
step functions of time in the data, and by forcing step functions by a straight- 
forward averaging of B 1 (t) and B 2 (t) to obtain < B t > and < B 2 > to be used in 
(I -1) is usually inadequate and possibly very inaccurate, as was discussed above. 

Part II discusses alternative methods, that is, multiple spacecraft methods, 
of obtaining accurate shock normals. It reviews established means by which our 
best estimate scheme can be tested, provided conditions are proper for the test. 
Part IV contains a short discussion of the use of the scheme in terms of simu- 
lated shock cases. And finally, the last section of Part IV deals with the actual 
calculation of a real shock normal previously studied by J. K. Chao (1970) and 
which serves as a test case of the overall scheme and associated computer pro- 
grams (which appear in Appendix C). 

It should be pointed out that this scheme accomplishes a good deal more than 
simply yield, in some sense, best estimate shock normals. It also provides best / 
estimate values for all the eleven relevant magnetic field and plasma parameters. 
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H. ALTERNATIVE METHODS 

By alternative methods of obtaining shock normals we refer to multiple 
spacecraft methods. These methods apply when the two or more spacecraft lie 
in a small enough region (R) of space such that the two following basic assump- 
tions apply: 

1. The shock surface can be assumed plane as the shock traverses R. 

2. The shock velocity remains constant over R. 

These limiting assumptions are, of course, not imposed by our single space- 
craft method described in Part XU below. We now briefly describe three multiple 
spacecraft methods of obtaining shock normals. 

A. Three Spacecraft Method 

Ogilvie and Burlaga (1969) and Greenstadt et al. (1970) employ the three 
spacecraft method to obtain shock normals. This method requires the shock to 
be observed at three spacecraft located in region R such that the plane in 
which the three spacecraft lie has a substantial shock normal component per- 
pendicular to it. Where R 12 and r 12 are the relative displacement vector and the 
shock time delay, respectively, between the first and second sightings etc., it 
is easy to show from simple geometrical arguments that, for the shock speed V s , 

r i 2 V s " ( n ' ^12) ’ 

r i3 V s = (n • R 13 ) , 
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and 


where. 



(n-i) 


such that 


n x 2 + n y 2 + n z 2 = 1 . (H - 2) . 

With these four equations the four unknowns n x , n y , n z and V s can easily be 
calculated. If the delay times are much longer than the uncertainties, in the 
time measurements, then this method is usually a very reliable one. 

B. Olbert Method 

Professor S. Olbert* of M. I. T. has devised a method that requires only two 
spacecraft observations, one of which needs only to record time of shock on- 
set and no other information (it could be the earth at sudden commencement). 

The other, however, must obtain magnetic field and plasma data as well as the 
shock onset time. By using the continuity of mass equation [See Equation (m-A~l)] 


•Private communication. 
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and the coplanarity theorem (Colburn and Sonett, 1966) Olbert shows that the 
shock normal is given by 

[b i x wj x U 

n " \fB ± x W) x U| (H- 3) 

where is the magnetic field before the shock, VI (= V 2 - V t )is the plasma 
velocity difference, and U is defined as 

R P 2^2 ~ p i ^1 

U " T P 2 - Pi ’ (H - 4) 

where R is the vector displacement and r the time delay between the two space- 
craft, and p x and P 2 are the plasma densities before and after the shock measured 
at one of the spacecraft (at which B 1} V 1} and V 2 are measured). 

This method is useful when reliable plasma data is available and when the 
magnetic field after the shock has relatively large fluctuations so that Equation 
(I - 1) can not be used. 

C. The Two Spacecraft Test 

If through some other method n is estimated, then the first of Equations 
(II - 1) constitutes a two spacecraft test of that estimate, provided V 3 can also 
be reliably calculated. This will not be a conclusive check but can serve as a 
means of "filtering out" some bad normal estimates and adding strength to the 
estimates of others. We will make use of this straightforward check in Section IY-C. 
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in. IMPROVED' SINGLE SPACECRAFT METHOD FOR 


OBLIQUE SHOCK NORMALS 

A. Theoretical Basis and Conservation Equations 

The shock normal improvement scheme described here rests on the follow- 
ing assumptions: 

1. The Rankin-Hugoniot conservation relations expressed for an isotropic 
medium are applicable to interplanetary type shocks (Ogilvie and 
Burlaga, 1969 and J. K. Chao, 1970). 

2. A shock can be represented as a "noised -up" step function increase in 
time as described in Section III-C, 

3. Magnetic field and plasma (proton) bulk velocity and density data provide 
adequate observational information for our purpose. That is, temperature, 
pressure, electron-data etc, are not necessary for significant normal 
improvement even though they might be necessary to strictly identify 

the shock in the first place. 

Only oblique shocks are considered in the scheme. That is, the special cases 
of the normal being either parallel with or perpendicular to the magnetic field 
are not treated here. 

We now begin by stating the basic equations of our system. 
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The conservation equations in the shock (*) frame of reference for an 


isotropic medium are 


H ]i ■ ° ■ • 

£> v „‘v; - = 0 

[ v »’B t - V.'Bn]’ = 0 , • 


where t » t x ort 2 



= 0 , 


P + (b 2 - B^jfzTT + pV* 2 2 


0 , 


and 


<m A - 1) 
(HI A - 2, 3) 

(HA - 4,5) 

1 

(m a - 6) 

pA-7) 


2 

= 0 , (III A - 8) 

i 

where p is the plasma mass density, V n * is the plasma bulk velocity component 
normal to the shock surface, V* (t = t x or t 2 ) are the components tangential 
to the shock surface, B n and B t ( t = 1 1 or t 2 ) are the associated normal and 
tangential components of the magnetic field, P is the total kinetic pressure, n is 


y *2 y P , B 2 ' (B • n) (V* • B) 
y~ + p + " 4 up (v* • n) 
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a unit vector normal to the shock surface, and y is the usual ratio of specific 
heats for the plasma. The symbol [ ] 2 means that the quantity within the 

brackets is to be evaluated before ("l") and after ( ,, 2 ,t ) the shock transition 
zone and then the quantities subtracted. Equation (HI .A - 1) is the mass con- 
tinuity equation, Equations (III A - 2 and 3) are the momentum conservation 
equations for the tangential components. Equations (in A - 4 and 5) are the 
tangential electric field continuity equations, Equation (El A - 6) is the normal 
magnetic field continuity equation, Equation (311 A - 7) is the momentum con- 
servation equation for the normal component, and finally Equation (m A - 8) is 
the energy conservation equation. According to assumption #3 above only the 
first six of these eight equations will be used in the normal improvement scheme. 
One sees that these equations can not be used directly without knowledge of 
n and the shock speed. Conversely then these equations may be viewed as 
constraints on the allowable values of n for a given set of relevant shock data. 

It is in this indirect sense that these equations will be used. 

B. Overdetermination Equations in Arbitrary Reference System 

The first six conservation equations, (III A - 1 to III A - 6), can be separated 
into two sets, three equations in each. We call these sets the shock velocity set 
and the overdetermination equations set. Appendix A demonstrates how this 
separation is made and provides a proof of the overdetermination equations. 

The shock velocity set is 

V s = V s n (3 equations) , (HI B - 1) 


11 




where 


{p 2 v 2 ~ Pi v i) ’ n 


P i P i 


(in b - 2 ) 


n s 


AB x (B t x B 2 ) 

(Bj x b 2 )| 


(m b - 3) 


AB = B 2 - B x , 


(m b - 4) 


and where the transformation equation 


V. = V; + V s (i = 1, 2) 


(in B - 5) 


was used. 

And the second set, constituting the remaining (three) overdetermination 
equations, is 



(HI B - 6) 


1 

1 

CM 


in • 

W x (b i X B 2 J 


= 0 , (IH B - 7) 
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and 


where 



0 (III B - 8) 


W a 



(HI B - 9) 


Firstly, we notice that Equations (HI B - 6, -7, -8) are rendered in general 
vector form and are independent of the shock (*) frame of reference. Therefore 
they can conveniently he used in association with whatever coordinate frame the 
experimenter wishes. Their simplest use then will be for a frame fixed to and 
moving with the measuring spacecraft and oriented in some physically meaning- 
ful way. The arbitrary system will have x - y - z axes by our terminology, 
where for instance the x axis might be along a direction radially away from the 
sun and the z axis normal to the ecliptic plane etc. According to this format any 
three dependent variables can be isolated through the use of the three equations. 
Chosing these to be p 1 ,p 2i and W x the overdetermination equations become 


"i 



5979.14 



R , 


(in B - 10) 


N 


2 



5979.14 (r - 1) R , 


(HE B - 11) 


13 



and 


W x = E , 


(III B - 12) 


where m p is the mass of the proton, and N 2 are in units of number of protons 
per cm 3 , all velocities are in km/see, and magnetic fields are in 7 , and where 


r 


/B„ S + B„ S + B„ S 

[ 2x x 2y y 2z z 

VB- S + B. S + B- S 

\ lx^x 1 y y 1 z z 


(HI B - 13) 


R = 


FG 
ID > 


E = 



r 

S s ff Q - W Q , 

x y ^z z v y ’ 

S = W 0 - E Q , 

y z v x- v z ’ 

S a E Q - W Q , 


(HI B - 14) 


(IHB- 15) 


(ELI B - 16) 


(El B - 17) 


(EE B - 18) 
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T = E x 0 + w y Y 0 + w z z 0 , (m B - 19) 

D h EM + V II. +1 M, , (HB-20) 

x y y 2 z 

G , + B 1J1 M y + B lz M z , (m B - 21) 

M, - Q,V 0 - Qy z 0 , (MB -22) 

M y . Q z Z 0 - Q z x 0 , (in B - 23) 

M « * <3, x o - Q, • (MB -24) 

Qx ■ B i s B 2* - B l z B 2 y • (mB-25) 

Qy - B I z B 2 x - B 1 x B 2 z . (IH B - 26) 

Q z = B ix B 2y ^ly®2x ’ (HI B - 27) 
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F 


(III B - 28) 


* ( x 0 2 + Y 0 2 * Z*)/«r . 


x o * 




(m B - 29) 


Y 


o 


B 0 - B- 

2y ly 


(III B - 30) 


and 


Z 


o 


B 2z " B 


1 z 


(HI B - 31) 

4? 


and also from Equation (III B - 9) 


W 


y 


V, - V, 

2y ly 


(m B - 32) 


and 


w 


3 V - V 

2z v lz 


(HI B - 33) 


That is, once eight parameters are fixed the remaining three, N t , N 2 and W x , 
are constrained to take the values dictated by Equations (III B - 10, -11, -12), 
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This constraint is the physical basis for the best-estimate scheme to be de- 
scribed in Section III - D. 

It should be pointed out here that if Equation (III A - 7) were written in the 
notation given above, it would become 


D 2 R A 

Ap ■ p 2 " p i = -W - W ' (HI B - S4) 

where the change in pressure across the shock surface, AP, is in units of 10“ 1 0 
dynes per cm 2 and where 


M 2 = M 2 + M 2 + M 2 

x y z 


(III B - 35) 


and 


A 


V* 


B 2y 


t> 2 _ n 2 „ 

B 2z 


■n 2 _ 

B iy 


B i 2 , 


(nr b - 36 ) 


Equation ( in B - 34) does not play a direct role in the estimation scheme but it 
can be used to calculate AP from the best-estimate parameters resulting from 
the scheme. Then the value of AP can be compared directly to pressure data 
(obviously the electron pressure cannot be ignored if this comparison is made). 

Since W = V 2 - V 1 , it is easy to see that Equation (in B - 2) can be written 
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as 

N 2 W • n 

v . = < nI B - 37 > 

where N i = p.J m p (i = 1, 2) was used. In Equation {HI B - 37) the first term 
and the factor n in the second term (from Equation (III B - 3)) can be readily 
calculated from the best-estimate results,' Then, in general, the calculation of 
the shock speed V s will be only as reliable as the value of Vj , the undisturbed 
pre-shock plasma velocity. However, using a straightforward average to obtain 
V x should give an adequate result, because the rms deviation on the magnitude 
of the pre-shock velocity is usually only a small fraction of the magnitude itself 
and its direction fluctuates very little (differing from the radial only slightly) . 
This depends somewhat on the provision that a proper averaging interval is 
chosen. Experience shows that a proper interval might characteristically be 
anywhere between 5 and 25 minutes. Finally, Equation (III B - 1) is used to 
obtain the vector shock velocity. 

C. The Noise Problem 

The usual conceptual model of an observation of an interplanetary type shock 
consists of a step function increase in time of the magnitudes of the shock quan- 1 
tities B, V,P (or T), and N as one goes from the upstream to the downstream 
positions. [For a so-called slow shock |b[ must decrease (J. K. Chao, 1970)]. 
The transition zone thickness is usually on the order of seconds, unless the probe 
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is observing the shock surface traveling edge-on. Indeed, for each physical 
entity the conservation relations accept only two values (a ’’before" and an 
"after" transition value). We retain this exceedingly simple concept of a shock 
but with the addition of stationary, uncorrelated, zero mean, noise to each of the 
basic shock quantities. That is, the noise is mathematically represented by a 
stationary, uncorrelated, zero-mean, random process. In most cases, however, 
we will find it necessary to restrict the before and after time zones to about 15 
and 10 minutes respectively. Other cases might require longer time zones. 
Figure 2 describes the shock model used in this work. Pressure (P) and tem- 
perature (T) are not shown because they are not used as part of the estimation 
scheme. 

MODEL OF SHOCK 



MINS. 


Figure 2. The conceptual model of an observation of an interplanetary type shock treated in this 

report. The straight line segments (i.e. the step function part) of the dashed curve refer 
to the basic (“true”), underlying vajues of the magnitude changes of the shock para- 
meters shown. For a slow shock |b] changes in the opposite direction to that shown. 
Time intervals are only approximate. 
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D. Least-Squares Loss-Function Procedure 


Guided by Section in - B we split the eleven basic parameters of our system 
into two groups: the 8 -parameter independent set and the 3 -parameter dependent 
set . The eight independent parameters then are chosen for convenience to be: 
B lx , B ly ,B lz ,B 2x , B 2yf B 2z , W y , andW z . Therefore, the three dependent 
parameters are N 1} N 2> and Vf x . The coordinate system x-y-z is an arbitrary 
orthogonal system and therefore in a R-T-N system (See Figure 1), for instance, 
W x may be W R or W T or W N etc. provided one is consistent with the use of B x and 

b 2 . ' 

In our mathematical scheme it will be convenient to define two vectors Y 
and X in the following fashion. The Y vector is the vector of observations, the 
so-called data array. If a total of N observations (including all data types: B lx , 
W y , etc) is to be used in the scheme, then Y will have dimension N. We have 
eleven basic data types and we impose on these types an order so that it becomes 
meaningful to speak of the first data type, the second data type, etc. Define 
N(i), i = 1, 2, . . . 11, as the number of observations of the i-th data type. 

Then the first N (1) elements of Y are to be the observations of the first data 
type, no particular order being necessary within the type, the elements of Y 
from N (1) + l to N(l) + N (2) are to be the observations of the second type, in 
any order within that type, etc. In symbolic form we write the N dimensional 
vector Y as 
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where, of course, 


1 1 

N = J2 N(i) • (HID -2) 

i=l 

Now we define the scalar symbol X i} a variable, to be the "best estimate" 
of the shock parameter measured by the observations in the i-th data type. The 
definition of what constitutes a best estimate will be provided below, and from this 
definition a mode of calculating a numerical value for X., i = 1, 2, ...11, 
will be evident. Define an N dimensional vector X by permitting each of the 
first N (1) elements of X to be Xj identically, the nextN (1) + 1 toN (1) + N(2) 
elements each to be X 2 identically, etc. Symbolically we write 


X 


« ( x i; Xy x 2 , ... X 2 , ... x n 




X tl ). (HID-3) 


repeated N(l) times N(2) times N(ll) times 


The X and Y arrays must be compatible component for component with respect 
to the parameter types, i»e., by the i-groupings. Notice for later reference that 
any component of Y can be expressed as 3 ^ 1 ^ where j ( i) refers to the j-th 



observation (or j-th point) of the i-th data type. Now for definiteness we make 
the following identification: 


Parameter Type 


r- 


Bi, , 
B, 


'ly 


B 


1 Z 


B 


2 x 


Table 1 

Data Array 
Component Symbol 

Y x < *<*>>, 3(1) = 1,2, ... 

y 2 (j(2)) 9 j(2 ) = 1, 2, ... 

Y 3 < j < 3 >\ j(3) = 1, 2, ... 

Y 4 (i < 4)) , 3(4) = 1, 2, ... 


N(l) 

N(2) 

N(3) 

N(4) 


Best estimate 
Array Component Symbol 

Xj repeated N(l) times 

X 2 repeated N(2) times 

X 3 repeated N(3) times 

X 4 repeated N(4) times 


independents 


r 

dependents 

v_ 


B 2z 

W y 

Wz 

N 1 

n 2 

w. 


Y s (j<S)) , 3(5) *1. 2, ... N(5) 
Y 6 < j < 6) \ 3(6) = 1, 2, ... N(6) 
Y 7 (i < 7)) , 3(7) *1, 2, ... N(7) 
Y 8 (j(8)) , 3(8) = 1. 2 » »• ^(8) 
Y 9 <j< 9 », j(9) = 1, 2, ... N(9) 
Y 10 ^‘< l °)>j(10) = 1, 2, ... N(10) 

3(11) = 1, 2, ... N(ll) 


X s repeated N(5) times 
X 6 repeated N(6) times 
X 7 repeated N(7) times 
X 8 repeated N(8) times 
X 9 repeated N(9) times 
X 10 repeated N(10) times 
X lx repeated N(ll) times 


This scheme will be used throughout the remainder of this work. 

The dependent parameters are related to the independent parameters through 
the overdetermination equations given by Equations (HI B -10, -11, -12). In the 
new notation these equations are formally expressed by the following: 


X 9 = X 9 (Z) , 
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X 


i o 


= x 10 (Z) 


X n = ' x u( Z ). 


(in d - 4) 


where 


Z - (x if x 3 , ... x 8 ). (HID - 5) 

That is, X 9 , X 10 , and X 1:l are functions of X 1( X 2 , ... X 8 only, rendering X, given by 
Equation (ED - 3); in terms of eight implicit variables, which must yet be 
determined. 

Also we define, the "sigma noise parameters" to be: 

cr. = a l + At (i = 1, 2, . . ‘ 11) (HI D - 6) 


where 


O -.' = 

1 


N(i) 


1/2 
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j=i 


N(i)-1 


(HI D - 7) 
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the unbiased rms deviation of the ith parameter-type observations, and 




N(i) 


r y . o! 
/ > 
j=i 

N(i) 


(HI D - 8) 


the average of the ith parameter-type observations, and where 

Act. (> 0) (HID -9) 

is additional weight given to cr.' to account for instrumental noise. 

We now define a sealer quantity known as a loss function, which is a 
measure of how well Z "fits" the data array Y. The smaller the loss function the 
better the fit. For this function we choose a standard cr -weighted least squares 
loss function: 


L(Z) 



(HI D - 10) 


Notice that L is a function of Z only, i.e. a function of only X lt X 2 , ... X g . Other 
functional forms for the loss function could be used provided they are positive 
definite. The exact structure of L is, of course, somewhat arbitrary. We define 
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the best estimate of Z to be the value of Z which minimizes the loss function 
<m D - 10). 

In order to minimize L(Z) its gradient with respect to X t , X 2 , ... X 8 must be 
zero. Hence, we set 


'~dXT = 0 ( i = 1 , 2, ... 8) (HID -11) 

for a necessary condition of solution. Because of the nonlinearity of Equations 
(HI D - 4), the eight equations given by Equations (m D - 11) represent a non- 
linear set to be solved simultaneously for the eight unknowns, the components 
of Z. Strictly speaking it is the solution of these equations which yields the 
components of the best estimate array. Expression (III D - 3) is more precisely 
a variable state vector whose all eleven components become, with the help of 
the overdetermination equations, the best fit array upon imposing condition 
(HID - 11). 

An iterative procedure will be used to solve the eight equation set (III D - 11). 

The numerical techique used is the Newton-Raphson method. See Appendix B.l 

♦ 2 ♦ ^ 

for a more detailed development of the overall statistical methods and the 
numerical technique in use here. We outline below the numerical procedure. 

We define Z as the exact solution of Equations (H[ D - 11) when an absolute 
minimum is attained. Then Z Q is defined as the first estimate (i.e. the "starting 
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vector" for the iteration procedure) of Z.* The vector Z 0 could be, for instance, 
the average of the first eight data channels of Y, i.e. its components could be 

<*!>.{*»>■ •••<*•>'• 

For AZ n . defined as 


AZ ■ Z - Z , 

n n n~ 1 ’ 


(in d - 12) 


Equation (B - 10) of Appendix B.l shows AZ n in explicit terms to be 


Az „ - [fiA T e;‘(Y-x)]|_^ 


(El D - 13) 


where 


A(Z) 


dX(Z) 

dz 


(HID.- 14) 


whose elements are A. . =5 dxJdX^ i = 1, 2, ...n[n is given by Equation 
(El D - 2)] and j = 1, 2, ... 8, • 


i(Z) - (a'q^aY 1 , 


(El D - 15) 


*In the strictest sense this should he Z defined in Appendix B. 1. But the statement is still 
correct, in a relative sense, as it stands* 


26 



and where . 



such that o-j 2 is repeated N (1) times, cr 2 2 repeated N(2) times, etc. By repeated 
application of Equation (HI D - 13) with 

K = Z„-l +S Z„ (DID -17) 


for n = 1, 2, 3, , provided Z Q is carefully chosen to insure convergence of L 

to its absolute minimum, Z n should tend toward Z, the exact solution of Equations 
(m D - 11). This iterative' procedure can be discontinued after a fixed number 
of steps or when | AzJ becomes sufficiently small, i. e. when 


|AZ„| 



< e 


(m D - 18) 


for some sufficiently small e > 0. We fix e at 0.01 and set n max , the total 
number of iterations allowed, equal to 15. The iterative procedure continues 
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from n = 1 through n = 15 unless criterion (IE D - 18) is satisfied. If a poor 
choice of Z 0 is made the process may diverge. In this case, except under very 
unusual circumstances, (El D - 18) will not be satisfied, and the process stops 
at n = 15. Then a new Z 0 must be chosen. Since the loss function can be cal- 
culated at each iteration step, then, even for a diverging case, that Z n associated 
with the smallest loss function is the one nearest to some acceptable starting 
vector Z 0 in a least squares sense. It must then be slightly changed in usually 
only a few components to provide an adequate Z 0 . Because of the nonlinearity 
of Equations (IE D - 4), and hence the nonlinearity of Equations (ni D - 11), the 
iterative process may converge to a false minimum, i.e. to one other than the 
absolute minimum sought. It is obvious when this occurs, because it leads to 
the "best estimate 11 values of Nj, N 2 , and W x differing greatly from the average 
values of these quantities, i. e. by more than 2 o- for one or more of the three. 
Other hints of a false convergence are results leading to N x > N 2 or > P 2 when 
the seventh conservation equation, in the form of Equation (El B - 34), is used. 
This false convergence also requires trying a new Z Q to bring about true 
convergence. 

In this connection it is useful to define a quality index, q, by the following 



(El D - 19) 
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where Z c is that value of Z which leads to convergence within some e, and N, 
given by Equation (in D - 2), is the total number of data points for all eleven 
parmeter-types. Obviously the nearer Z c is to Z the larger q will be. The 
quantity q should be near, or slightly greater than, unity for common cases of 
interplanetary shocks. For too small a q, say q = 1/2 or so, the convergence 
may be a false one , 
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IV. TEST OF METHOD AND EXAMPLES OF ITS USE 


A, Error Cones and Simulation 

The preceding section outlined a method of utilizing plasma data to obtain 
"good" estimates of the before and after magnetic fields. Another and far simpler 
method of estimating these fields is to take straightforward averages of the ob- 
servations of the fields as the estimate. The justification for utilizing the more 
complicated weighted least squares estimation procedure instead of the mean 
value method has been that the more complicated method yields a more accurate 
estimate of the before and after fields in general. And any small increase in the 
accuracy of these estimates because of the form of Equation (El B - 3) can yield 
substantial improvements in the estimate of the shock normal. 

But, of course, this is an assumption which must be tested and proved, at 
least within some reasonable basic set of assumptions. In short, it is necessary 
to show that the weighted least squares estimation procedure leads to significantly 
better estimates of the shock normal than the mean value proceedure within the 
limiting assumptions stated in Section m - A. In Section IV-B we do this by 
applying both estimation methods to simulated observations of a shock and 
associating with each method ah error cone about the true normal to the shock 
surface. The comparison of these error cones will indicate the degree of im- 
provement to be expected from the weighted least squares technique. 

The meaning of these error cones and the means by which they are calculated 
will now be described. First, simulated shock observations are generated by 
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assuming a simple underlying step' function model of 'a shock (See Section Ift-C) 
and associating with such a shock the eight independent parameters which con- 
stitute the components of the vector Z. These parameters are chosen to be con- 
sistent witfi previous studies of shock properties. Then the remaining three • 

■ dependent parameters are obtained from the overdetermination equations of ■' - 
Section m-B. These eleven parameters constitute the "true" shock parameters 
of the simulated shock. Zero mean (i.e. unbiased), stationary, uncorrelated, 
normally distributed noise is then imposed on-all of the eleven measurable 
parameters. The number of observations of each measurable parameter and the 
variance- figures on the noise-are again'chosen to be typical of what one should 
expect of shock observations. ‘ Both the weighted least squares and the mean 
value techniques are then utilized to obtain estimates of the before and after - 
magnetic fields (and estimates of the plasma parameters). It is possible-to 
obtain covariance -matricies for both these estimation procedures. The manner 
by which this is done for the least squares method, along with general mathe- 
matical details of error cone construction, is given in Appendix B.2. The . 
covariance matrix of the mean value estimate can be easily obtained by .recalling 
that-the variance of an estimate obtained by a mean value- is just the variance of 
the underlying population divided by the. sample -size. This provides us with -the 
diagonal elements of the desired covariance matrix. And since the noise on 
each data type -is assumed to-be independent of the noise on other data types, 
the off diagonal elements are zero. 
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These matrices are measures -of the statistical dispersions of the estimates 
about the true values. "What interests us now is how these statistical dispersions 
propagate their way through the non-linear function, of Expression (IE B -.3) into 
angular errors in the estimates of the true shock normal! Specifically we shall 
obtain, for each estimation procedure, a 95% critical angular error value a, that 
is, an angle for which the probability of the angular error (caused by the use of a 
particular estimation procedure) being smaller than a is 0.95. A Monte Carlo 
process is necessary .to obtain these critical angles. Essentially this Monte Carlo 
procedure represents a method, indeed the only method, of propagating the statis- 
tical dispersion of the magnetic field estimates, as measured by a covariance 
matrix, through the highly non-linear function (IE B - 3). The resulting critical 
error angle is related to the statistical dispersion of the. estimate of the shock 
normal. The Monte Carlo procedure is also described in Appendix B.2.- 

These critical angular error values have an obvious geometric interpretation, 
namely, a 95% critical error angle a can be represented by the defining angle of 
a right circular cone with its axis being the true shock normal. 

In Section IV-C error cone angles associated with a real shock will be 
calculated. True normals are not available in the cases of real shocks, of course. 
Hence, in these cases the two methods of obtaining error cones, the mean value 
and least squares methods, must have error cones defined in a slightly modified 
way from those of the simulated cases. In the mean value case the axis of the 
error cone will be the normal obtained from taking straightforward averages of 
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the field, and the cone is generated with respect to this average normal.. In the 
best estimate (least squares) case the axis of the error cone is the best estimate 
normal, and the cone is generated with respect to that normal. It must be 
emphasized that these real shock associated- error cones can not enjoy the same 
rigorous interpretation as those of the simulated shocks. But since the cone 
angles are expected to depend strongly on the cr. ’ s (i = 1, 2 ... 11), defined by 
Equation (IQ D - 6), and only weakly on the actual shock parameters in most 
cases, then the real cones, for practical purposes, should have interpretations 
analogous to the cones of the simulated cases. That is, the probability of the 
true shock normal lying within the 95% critical error cone for real cases is 
approximately 0.95., As in simulated cases of shocks the best estimate error 
cone for a real case will have a cone angle smaller, and sometimes very sub- 
stantially smaller, than the mean value cone angle. 

B. Study of a Simulated Case 

As described in Section IV-A and Appendix B.2 realistic simulated shocks 
were generated in order to test the degree of success of the improvement scheme 
and to provide a check on the associated computer programs. The program has 
the capability (See the XMONTE subroutine in Appendix C.l) of generating a 
simulated Y data array using preassigned values of a i , N(i) (i = 1, 2, ... 11), and 
X L (i = 1, 2, ... 8), the latter being components of what we refer to as the Z true 
vector*. The ’’true" components X ».X 10 , and Xj j (dependent parameters) are 

*This is called Z in Appendix B. 
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obtained from the overdetermination equations of Section m-B. The eleven 
true X i ’s are "noised-up" by adding the output from a random number generator, 
according the values of the cr. *s and the N(i)*s, yielding the vector Y . Con- 
sistent with our previously described model the random number generator 
provides, for all practical purposes, samples of an unbiased, stationary, un- 
correlated, normally distributed random variable. By calculating the mean 
values of the first eight data types, using Y as if it were real data, gives X^mean) 
i = 1, 2, ... 8, or Z M in vector notation. The Z M vector should, in most cases, 
provide an acceptable starting vector for the iteration procedure of Section ni-D. 
We then set Z Q = Z M for all simulated shocks. Hence, we enter the simulation 
problem with all of the information that would be necessary to employ the im- 
provement scheme to a real shock with the important difference that here Z true 
is known. And by design, the simulated data does satisfy the statistical model. 

Kealistic input parameters cr it N(i), (i = 1, 2, ... 11) and Z true were used 
to test the program. Table 2 gives an example of input values used in such a 
test. The shock computer program is listed in Appendix C.l. Appendix D.l 
shows an example of a printed output of the results of using the input values 
given in Table 2. It represents only one 140-number sample from the random 
number generator, where 140 is the sum of the N(i) T s. Any number of samples 
from the generator, each giving a different Y , are available where, of course, 
each Y represents just a single data sample of the true shock of Table 2. The 
preface of Appendix C.l explains what switches have to be set, and to what 
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Table 2 


Example of Input Values for Simulated Shock Test. 


Parameter 

i 

N(i) 

cr. (in units 
of X x ) 


Xj (true) 


i 

20 

0.35 

4.0 

*\ 

7 


Bly 

2 

20 

0.50 

5.0 

y 


El* 

3 

20 

0.35 

-1.0 

y 


B 2x 

4 

10 

0.60 

3.5 

y 

^independent 

B 2 y 

5 

10- 

1.10 

9.0 

y 


B 2 z 

6 

10 

1.30 

-3.0 

y 


W y 

7 

10 

10.0 

10.0 

km/sec 


w z 

8 

10 

10.0 

20.0 

km/ sec 


Nx 

9 

10 

0.7 

(7.26 

f/cm 3 ) 


n 2 

10 

10 

1.0 

(13.86 #/cm 3 ) 

^-dependent* 

w x 

11 

10 

10.0 

(75.83 km/see) 



•Strictly speaking the three values in. parenthesis are not input parameters. 

values, in order to run a simulated shock program (and also for a real shock 
program). In the particular case of this- 'simulated example the switches were 
set to the following values: 

IPRp = ' 1 
ISWTCH = 1 

ICASE = 2 

ISAMPL = 5 
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IPRO equal to 1 means simulated shocks being processed, and ISWTCH equal 
to 1 means that XSTART (which is the same as input Z true for a simulated 
shock) is replaced by XMEAN as the starting vector for the iteration process. 
ICASE of 2 simply means that two basic input shocks (or two cases) are being 
studied, where here Table 2 gives the values for only one of the two cases. 

ISAMPL equal to 5 means that five samples of noise are to be imposed in 
each of the two cases creating ten Y's. The input arrangement is such that the 
same basic shock is associated with each of the five samples, for both cases; 
that is, for each of the first five noise samples the same values of Table 2, say, 
are used as input to the random number generator and for the second five samples 
the same values of some other table (not shown) are used as input to the generator. 
Our sample output, Appendix D.l, is then the result of one -of the ten Y f s. Below 
we describe the shock program output sample. 

From what has been said above the first eight lines are self evident (where SIG 
is o\ and NNis N(i)). INPUT XSTART is justZ true . The so-called G values are 
the values of the quantities available from the overdetermination equations. 
Equation (HI B - 3), and Equation (HI B - 34) all of Section III-B (See the CON 
subroutine in Appendix C.l). These are the three dependent shock parameters 
Nj, N 2 , and W x , the x-, y-, and z- components of n, and the total kinetic pressure 
change AP. "Corresponding G values" then refer to those G values corresponding 
to Z true . The best estimate independent parameter matrix is a two dimensional 
array whose columns are the eight independent shock parameters, the components 
of Z, and whose rows correspond to the iteration steps. The top, or T, M - row" 
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(see far right for label), is composed of the mean values of Z, and immediately 
below that is the starting vector Z 0 in the zeroth-row. Notice the Mth and zeroth 
rows have the same values because XSTAJRT is replaced by XMEAN (which is 
not the case for a real shock). The process went the full 15 steps because 
Z = |AZJ/| Z n _ J did not become < 0.01 as Expression (m D - 18) requires for 
a number of iterations (L) less than 15. The Z ratio is printed out at the far 
right, and the value of the loss function also, at each iteration step. Below that 
the quality, defined by Equation (m D - 19), is printed out for each step. Below 
the independent parameter matrix is the associated dependent parameter matrix 
whose columns are the values of the G-parameters described above; the rows 
again correspond to the iteration steps. The B matrix is the evaluation of B , 
given by Equation (in D - 15), for Z = Z, shown as the last step of the independent 
parameter matrix (i.e. the best estimate step). The contracted derivative matrix 
A Is a contraction of A, given by Equation (HI D - 14), evaluated for Z = Z. By 
contracting A no information is essentially lost (See the C matrix of the AA 
subroutine in Appendix C.l); the statement just below Equation (in D - 14) 
concerning the elements of A, along with Equation (IE D - 3), explains why this 
is so. The three numbers at the very bottom of the printed output refer to angles 
in degrees. These are: 

AAVE is the angle between n (true) and n (Mean), 

ABE is the angle between n (true) andn (B. E„), 
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and 


AVE, ABE is the angle between n (Mean) andn(B. E.). 

Note: In the real shock output n (true) is replaced by n calculated from 
Z 0 , i. e., from XSTART. Hence, of the three angles only AVE, ABE has 
any relevant meaning in a real case , 

The top or M-th row of the independent parameter matrix contains all eleven 
mean values of the shock parameters, the last three of which are, strictly speak- 
ing, not part of the matrix. The first eight parameters are the ones used to 
obtain what below are called the MEAN’S G’s, which are self explanatory. Notice 
that the mean values N 1} N 2 , and W x are distinctly different from those same 
quantities derived from the MEAN'S G’s, and this is most important (especially 
in real shock cases). When this difference is very great it indicates the low 
quality of using straightfoward mean values as final estimates for the shock 
parameters. In fact, the mean quality parameter, QUALITY M = 0.127, is quite 
low compared to unity or so, which is expected for a least squares best estimate. 

Notice that after only about 4 or 5 steps the calculation is essentially com- 
pleted, and little is gained after those steps. The choice of e = 0.01 in expression 
(HI D - 18) is obviously a conservative one since this sample output is rather 
typical. A comparison of the true shock parameters (i. e. XSTART and cor- 
responding G’s) with those from the mean value and best estimate calculations 
shows for this case, or rather for this sample of a case, how valuable the scheme 
can be. But the true test of improvement lies in a comparison of the two methods 
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of obtaining error cones as described in Section IV-A; below this is discussed. 
Notice that the best estimate normal lies only 3 6 away from n(true) but n(Mean) 
is almost 9° away. In some cases (i.e. for other X^trueJ’s, etc.) the improve- 
ment might be much more dramatic and yet in others the improvement is in- 
significant. It is even possible that ABE turn out to be larger than AAVE, as one 
should expect in a statistical problem of this type, but it must occur infrequently. 

Appendix D.2 shows a sample of the printed output of the cone computer 
program listed in Appendix C.2. This example corresponds exactly to the shock 
case described by Table 2, and, of course, gives the cone angles associated with 
the example program output of Appendix D.l (and all other samples of this same 
case). The so-called "FI 1 RESULTS" refers to cone angles found by using the 
least squares technique, and the "FI 2 RESULTS" refers to the mean value method. 
All angles are given in degrees. The designations 15-, 30-, and 150- VALUE 
refer to 99.5%, 99%, and 95% error cone angles, respectively, where a Monte - 
Carlo sample size of 3,000 was used. For example, consider the 150 VALUE 
case: 3,000 - 150 (= 2,850) refers to 95% of 3,000, and designates the cone with- 
in which 95% of the normal estimates lie. We will not be concerned with the 
99.5% and the 99% error cones in this study. Notice then that, in this case, the ' 
angles a(Mean) and a (Best Estimate) are 10.6° and 5.2°, respectively. This 
represents an error cone angle improvement of better than a factor of 2, and is 
characteristic of realistic cases in general or perhaps is somewhat conservative. 
Sometimes the improvement factor is more dramatic (i.e. values of 3 and 4) for 
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realistically simulated shocks, and in no realistic case will it turn out that 
a (Mean) £ a (Best Estimate). 

C. Example of Actual Case; The August 29, 1966 (Pioneer 7) Shock 

The August 29, 1966 shock, as observed in Pioneer 7 data, was first studied 
by J. K. Chao (1970) and is reexamined here as an example of the use of the 
least-square technique described in this report. Taylor (1968) also observed 
this shock in the magnetic field data of Explorer 28 , but the associated plasma 
data was not existent for that spacecraft. 

In applying the least-squares scheme, 25 alternate data points, representing 
12.5 minutes, were used for the B t field in the Y array of Equation (HI D - 1), and 

i 

\ 

18 points, representing 9 minutes, for the B 2 field, and 5 points, representing 
**8 minutes, for each plasma parameter, before and after the shock, were used. 
Thp quality index for the best estimate convergence value, as defined by Equation 
(HI D - 19), was 1.03, which is a common sort of value for interplanetary type 
shocks. For the total 154 data points this corresponds to a loss function value 
of 145. Table 3 gives the values of the shock’s relevant parameters, as well as 
the observed onset time. Thecr's for the magnetic field were obtained directly 
from a calculation of the rms deviation in the data. The cr’s for the plasma param- 
eters were found likewise with the addition of instrumental A o-’s [See Equation 
(IE D- 6)] to the statistical values. The values for the components of W(= V 2 - V 1 ) 
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Table 3 


Pioneer 7 Shock Event of 
August 29, 1966 (14:16:57.4±.8 U.T.) 





Best 

Parameter 

Average 

a 

Estimate 


Value 


Value 

B « oo 

-2.3 

0.57 

-2.30 

B 1X 

0.9 

0.65 

0.68 

B in 

-2.3 

0.35 

-2.27 

®2R 

-3*9 

0.70 

-3.70 


2.1 

1.7 

2.89 

B 2N 

-6.8 

1.5 

-6.98 

W R (km/sec) 

79.4 

6.90 

78.9 

W T 

25.2 

10.2 

27.9 

W 

-12.9(-166) 

7.40 

-17.8 

Nj (I/cm 3 ) 

4.6(0.098) 

0.46 

4.88 

n 2 

14.9(0.206) 

1.80 

13.6 

n R 

0.94 


0.945 

n t 

-0.06 


0.296 

n N 

-0.35 


-0.142 

AP (iff-.. SE2t\ 

\ CHI 2 / 

2.2 (-1.6) 

0.5? 

6.9 

Error Cone Angle 

25.3° 


6.0° 
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were obtained by performing a mirror image subtraction about the shock transi- 
tion time, of V t from V 2 , rather than a chronological one, as shown below: 


10 9 8 7 6 


Shock 


Time order of points 


Jump 


1 


5 4 3 2 1 

H 1 1 i h 


Mirror Subtraction 
\fj( 1) = y($) _ y (51 
W<2) = y(7) _ y(4) 


Chronological Subtraction 
wO) = yOO) _ y(5) 
y/(2) = y(9) _ y(4) 


ft<5) =y(10) Jy{1) 


^(5) n y( 6 ) — y(l) 


This yielded the smallest cr's for W. [This variation of cr with the choice of the 
manner of subtracting Vj from V 2 represents a slight violation of the ideal step 
model of the shock]. 

The average values taken directly from the data are given in the Table for a 
comparison with the best estimate results. The quantities in parentheses in the 
average value column are the values one obtains by using the average values of 
the eight independent parameters with the overdetermination equations of Section 
IH-B. Notice that the average values of W K , N x , and N 2 correspond poorly 
with those values calculated via the overdetermination equations. This is 
true even though the best estimate values, which satisfy the overdetermination 
equations exactly, and the average values do not differ very appreciably 
except perhaps in the case of the B 2T parameter. This demonstrates the sensitive 
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nature of these equations. In a similar manner Equation (HI B - 34) is used 
to calculate AP. When average value parameters are used Ap is seen to be . 
negative which is impossible for an actual shock. The average value of 2.2 x 10" 10 
dynes/cm 2 corresponds to the change in proton pressure only but the best estimate 
value of 6.9 x 10” 1 °- dynes/cm 2 refers to all particle species including, of course, 
electrons, and is expectedly larger. The angle between the n-average and n-best 
estimate is about 24°. Since great confidence is placed on the best estimate value 
owing to the large error cone angle decrease (25.3° to 6.0°) of a factor of about 4, 
then n calculated via average magnetic field values only would have given an 
unacceptable result. Also J. K. Chao (1970) finds a value for n(n R = 0.97, 
n T = 0.25, n N = -0.04) which differs from our best estimate value by less than 
7°. He also uses a best-fit technique (of trial-and-error fitting to the conserva- 
tion equations) to obtain the normal. 

Further evidence that n-best estimate is a dramatic improvement over 
n-average, in this case, lies in applying the two-spacecraft test described in 
Section II - C. This was done by utilizing the shock onset information obtained 
from Explorer 33 (See Figure 3), which also observed the August 29 shock. In 
an R-T-N' coordinate system centered at the earth the position coordinates of 
Pioneer 7 and Explorer 33 were, respectively, 

R 7 = (257, 119, 7.7) R e , 
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Figure '3, Explorer 33 and Pioneer 7 projected (onto ecliptic plane) positions during their observations of the August 29, 1966 shock. 
Theoretical locations of the unaberrated earth’s bow shock and tail are shown, as well as the edge-on view of the local 
plane surface of the August shock with normal rf. R is the relative position vector between the spacecraft. 



and 


R 33 = (54.9, 26.4, - 17.0) R E , 

where S £ is the earth's radius. This yields a relative position displacement 
vector R = (1.29, 0.59, 0.16) in units of 10 6 km. Since the onset time at Explorer 
33 was 13:28.5±.7 the delay time between sightings, t, was 48.5 minutes. Using 
the first of Equations (II - 1), where the trial n is our best estimate value, n, we 
obtain an "observed" V s , which is 


n • R 


s, obs 


= 471 


km 

sec 


For a "calculated" V s we use Equation (IH B - 37), where Vj will be simply 
the pre-shock average velocity and all other quantities best estimate ones from 
Table 3. This yields 


s , c a 1 c 


= 467 


km 

sec 


where 


Vj = (353, 13.7, 24.8) — • 

We see that the observed and calculated values of V s differ by less than 1%. 
This fine correspondence is partly fortuitous since the second term in Equation 
(HI B - 37) ( Vj * n), the weak link in the argument, is probably in error by 
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slightly more than 1 %, However, we believe that in this example, and in any 
ease similar to it, the least squares method of calculating the shock normal 
leads to a significantly improved estimate of the normal as well. as of the 
eleven relevant shock parameters. 
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APPENDIX A 


DERIVATION OF THE OVERDETERMINATION EQUATIONS 


Let II denote the plane containing B x and B 2 and define the unit vectors 


AB 

I All 


(A -1) 


t 


2 


B 1 xB 2 
ilj xb 2 


(A -2) 


and 


n =' tjXtj 


(A -3) 


where 


AB s B 2 - Bj . 


(A -4) 


Since 


t 


1 


= 0, 


{A -5) 
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then 


ABx (Bj x B 2 ) 

n |^x(b iX b 2 )|- ( A “ 6 ) 

We see that t 2 is normal to the n plane: Let 2 denote a plane perpendicular to 
II and containing AB . Hence, both 1 1 and t 2 must lie in the 2 plane. Then n 
(sTjXtj,) is a unit vector normal to the 2 plane, provided neither tj nor t 2 is 
zero. It follows that n lies in the E plane. The coplanarity theorem (Colburn 
and Sonett, 1966) demands that the shock plane's normal lie in the n plane. 

£ Notice that AB (or t j) is common to both the n and 2 planesj. One sees im- 
mediately, according to these definitions, that B 1 • n and ,B 2 • n are equal, as 
required by Equation (m A - 6) if n represents the shock surface normal. We 
are then justified in uniquely identifying 2 with the shock surface, t 2 and t 2 as 
tangential to it while perpendicular to each other, and n as normal to it. The 
situation is shown in Figure 4. 

Using Equations (A - 1 to 6) we can rewrite Equations (HE A - 1 to 6), which 
become 


{P\ v i* 'Pi v 2*) * n = 0 ' ’ 


(A -7) 



(A -8) 
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Figure 4. The II and 2(shock) planes perpendicular to each other and intersecting along theT-j direction. The vectors B 2 , and n 
lie in the II plane, andT^ and lie in the 2 plane, with iTnormal tojt. The 2plane separates the undisturbed pre-shock 
medium (“l”) from the post-shock medium ( 4< 2 ,f ). Note that AB= B 2 “ B y 


(A -9) 



V 1 *xB 1 -V 2 *xB 2 


ab = 0 


(A - 10) 


(v 


x B[ ' 


V 2 xB 2i 


(BixB 2 ) = 


(A - 11) 


and 


AB • n = 0 , 


(A - 12) 


where 


AV* = V * - V * , 


(A - 13) 


and where Equations (A - 7 and - 12) aided in obtaining Equations (A - 8 and - 9), 
We define a new velocity V related to V* , the plasma velocity as measured in the 
shock frame of reference, by 


V i = V + V s n 


(i = 1, 2) 


(A - 14) 


where V s is the speed of the shock frame, fixed to the shock surface, measured 
with respect to whatever frame V is measured in (which could be the spacecraft 
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the equation becomes 


W • (b i x B 2 ) = 0 (A - 17) 

where Equation (A - 14) was used. 

Now consider Equation (A - 9), which is, upon expansion, 

Pl {V •")[*• (b. *«,)] - \r [ss • (b, * B,)] = 0 

where (A - 16) was used. From Equation (A - 17) this immediately reduces to 
Ab • ( B 1 x B 2 ) = 0 since Bj • n ^ 0 in the cases that we are considering. But 

this is already expressed by Equation (A - 5) and therefore reduces to another 
identity of no further use to us here. 

Consider Equation (A - 11) now. By using Equation (A - 14) it becomes 

VVVV V s("^)] • (b!*b 2 ) = 0 . (A -18) 

The third term in the brackets, with the help of Equation (A - 15), is 

n x AB . 


\ P 2 ^2 “ P\ ^l) n 
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frame of reference). By using Equation (A - 14) together with Equation (A - 7) 
V s can be shown to be 


(/ ? 2 ^2 P\ ^i) * n 


(A - 15) 


% V 

We will now use Equations (A - 14 and - 15), with n defined by either Equation 
(A - 3) or (A - 6), to render Equations (A - 8 to - 12) in terms of V instead of V*. ' 
Clearly Equation (A - 12) represents an identity when n is replaced by the 
Expression (A - 6). In this sense it is not an ” overdetermined equation 1 ' and can 
not be retained as such. Next we consider Equation (A - 10) which is, after 
expansion, 

Yj* x B t ■ B 2 - V; x B, • B x - V 2 * x B 2 • B 2 + V,* x B, •' B t = 0 

or, by the operation exchange rule for triple scaler products (op rule), is 

v; * Bj x B 2 - V* ■ BjXBj- v 2 * ♦ B 2 X B 2 + V 2 * * B 2 X Bj = 0 . 

Noticing that the second and third terms are zero and defining 

W - V 2 - V, (A - 16) 
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Using Equations (A - 1), (A - 2), and {A - 3) and the op rule this becomes 


(/^2 ^2 ^2 X ® 1 Vj X B 


>) ■ i >] 


p 2 ~ Pi 


Replacing this back into Equation (A - 18) and noticing that t 2 • ( B x x B 2 ) 
- | B x x B 2 | ; we finally obtain 




0 


where r is defined as 


r 


Pi 

p\ 


(A - 19) 


(A - 20) 


and. where the op mile was again employed. 

Only Equation (A - 8) remains to be reduced. By the usual substitutions it 
can be written 


P 



h-V s ) (w ab)- 


B, 


4 TT 


AB 1 


0 • 


Replacing V s and n by Equations (A - 15) and (A - 8) respectively this becomes 


P_iPi 
P 2 ~ P\ 


W(W • AB) 



0 . (A - 21) 
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We refer to Equations (A - 17, - 19, and - 21) as the Overdetermination 
Equations because, in the sense that all eleven shock parameters of the system 
are assumed measured, any three parameters are overdetermined by these . 

. equations using the other eight (independent) parameters. Notice that these 
equations do not depend on finding the directions tj, t 2 , or n or on any para- 
meter depending on the shock (*) frame of reference, such as V*, as Equations 
(HE A - 1 to - 6) did, and even as Equations (A - 7 to - 12) did in part. 

In review then we see that Equation (A - 7) provided V s , Equation (A - 12) 
and Equation (A - 9) (through the coplanarity theorem) gave us the direction n , 
and the remaining three Equations (A - 8, - 10, and - 11), properly transformed, 
yield the Overdetermination Equations. 
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APPENDIX B 


THE LEAST SQUARES ESTIMATOR OF A 
SHOCK NORMAL AND ITS ASSOCIATED ERROR CONE ANGLE 

B.l The Estimation Procedure 

Below we develop in somewhat general terms the estimation procedure, 
necessary for a better understanding' of the less statistically oriented Section 
m-D. The dimensionality of all vectors and matrices is evident from the dis- 
cussions in Section III-D, and the notation used here is consistent with that 
Section. 

*\j 

Let the vector Z represent a state which is to be estimated.* Its components 
are to be conceptually identified with some (i, e. any eight) of the magnetic field 
and plasma quantities shown in Figure 2 as the shock’s underlying step function 
(denoted by dashed straight lines) and discussed in Section III-C. Let the vector 

r\j 'X* 

Y, of higher dimensionality (N) than Z, be another state functionally connected to 
Zby a known function X. Thus, Y = X(Z) . Assume that the vector v represents 
-a multivariate normal distribution with mean zero and a known covariance matrix 
0 defined by 

Q y s cpy(Y) '= e[(y-E(Y>) (y-E(Y)) T ]‘ - (B - 1) 

where E is the expectation operator [E(Y) then simply being the mean value of 
Y, <Y>] , and where the superscript T represents the transpose of the vector. 

*Z represents Z true . 
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The problem to be solved is the following: given one sample Y from the dis- 

'Vi _ 

tribution Y + v obtain, in some sense to be discussed below, a '’best" estimate of 
Z. In practice Y is considered a state which is directly observable and which 
has a known relationship to a state Z, the quantity to be estimated. The random 

_ 'V 

variable v should be thought of as the noise on the observation of Y, caused, by 
instrumental inaccuracies and natural but unexplained fluctuations in the values 
of the relevant parameters, i, e. unexplained in terms of the "known relationship" 
mentioned above. The assumption that v has normally distributed components 
with zero means is a valuable convenience from a mathematical point of view. 

But it has more to recomend it than mathematical convenience. Giving v a zero 
mean implies we have assumed that all systematic or modeling errors have been 
removed from the analysis. If significant modeling errors have not been removed, 
then no estimation proceedure is likely to provide an acceptable estimate of Z. 
Hence, little appears to be lost in assuming a zero mean for v. The justification 
for modeling the noise v as a normal random variable rests on the vague meta- 
statistical analogue to the law of large numbers which can be stated as follows : 

"If a large number of random variables are combined in a reasonably complicated 
fashion to form a single random variable, then it is likely that this random vari- 
able will have a nearly normal distribution." The assumptions of this m eta- 
statistical principal are usually satisfied when one is making observations in 
nature. Thus, the assumption that v is normally distributed has at least some 
reasonable support. 
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It remains to be clarified in just what "best” sense Z.is to be estimated. 
One such common estimation procedure insists that an estimator of Z be 
chosen such that the weighted sum of the squares of the differences between 
observed and expected observations is minimized. More specifically we define 
the so-called loss function L as 

L(Z) = (X(Z)-y) t ^(X(Z)~ y) 

where the vector Z is an independent variable which tends toward Z, defined as the 
best estimate of Z, as L tends toward an absolute minimum, and W is the so- 
called weighting matrix of the loss function. is generally set equal to Q ~ 1 or 
some slight modification of it. Then for our purposes the loss function, which is 
given as Equation (III D - 10), is 


L = (X-Y) t 2; 1 (X-Y) 


(B-2) 


where £> y is given formally by definition (B - 1). 

We minimize L in the following way, known as the Newton-Baphson method: 
The gradient, G, of Equation (B - 2) is 


G(Z) 


<9L 

a 



(X-Y) . 


(B-3) 
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We define a matrix A as the following 


A(Z) 


dX(Z) 

dZ 


(B-4) 


[Note that the elements of £ are . } ~ dX^dX^ , i = 1, 2, ... N (n given by 
Equation (m D - 2) }• and j - 1, 2, ... 8.] 


Prom Equations (B - 3) and (B - 4) the gradient of G is 


5G 

dZ 


d£ J 

2 - ^=‘ S y ~ 1 (X ~ Y) + 2£ r Q^£ . 


(B-5) 


The first order Taylor’s expansion of G(Z) is 


G 




t <?G(Z) 
+ — 

dZ 





(B-6) 


Using Equations (B - 3) and (B - 5) and disregarding the latter’s second order 
term Equation (B - 6) becomes 



= 2A T Q ~ 1 (X- Y)|-£ + 2A t O"1a|- AZ , 

_ n- 1 — — y n~ 1 


(B-7) 
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where 



Z 


n 


- Z n-1 


(B-8) 


For minimization of L(Z-), G^Z n j must be zero, as expressed by Equation 
(m D - 11). Then by defining the eight by eight matrix B in the following way 


B 



(B-9) 


Equation (B - 7) , upon minimization, becomes 


(B - 10) 


Combining Equations (B - 8) and (B - 10) yields 



- Z + I(Z)A T (Z) Q' 1 (y - X(Z)) 




(n = 1, 2, ...) . (B - 11) 


By repeated application of Equation (B - 11), provided Z Q is carefully chosen to 
insure convergence to the absolute minimum of L, Z n should tend toward Z, the 
exact solution of Equation (IQ D - 11). This iterative process should converge 

A ^ 

rapidly to the correct value Z, the best estimate ofZ . 
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For different samples Y of the random variable Y + v the iteration pro- 

A 

cedure will yield different values of Z, representing samples of the associated 
distribution of a new random variable. This new entity will henceforth also be 
symbolically represented as Z since there is little chance iof confusing the 

o 

random variable with one of its samples, i. e., the solution of G(Z) = 0 for a 
given Y . 

A *Vr 

In evaluating the quality of an estimator Z of Z , two factors are usually 
considered. One is the extent to which the estimator is biased. The bias of the 
estimator Z of Z is defined as E[Z - Z] . The other factor is the dispersion of 

A f\s A 

the estimator Z about Z . This is obtained by taking the second moment of Z aboul 

/\j f*" A A r\f j 

Z defined as E (Z - Z) (Z - Z) T ! , Of course, the smaller the bias and the dis- 
persion the better the quality of the estimator. Neither the bias nor the dis- 

A 

persion of Z can be conveniently calculated without the imposition of a certain 
linearity assumption. It must be assumed that X can be represented by a linear 
expansion of itself about Z . From Equation (B - 4) this is 


X(Z) - X(Z) = A&) (Z - Z) . (B - 12) 


If one assumes that the vector root Z of G(Z) = 0 is sufficiently close to Z to 

A 

permit the use of Equation (B - 12), then G(Z) = 0 can be written as 

A T (Z)fi y _1 X(Z) + A T GOQ^ACZ) (Z - Z) - A r (Z)^" 1 Y = 0 , (B - 13) 
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with the help of Equation (B - 8) . By applying the expectation operator to ' 
Equation (B - 13) one obtains 


(Z) 2 y " 1 X(Z) + A X (Z)Q ~ 1 A(Z ) E[Z - Z] - | T (Z) £ y “ 1 E[Y] = 0 . (B - 14) 


Since Y is a sample chosen from the distribution Y + v = X(Z) + v and since 
E[v] = 0, as mentioned above, then it follows that E[Y] = X(Z). Therefore 
Equation (B - 14) yields : 


E[Z - Z] = 0 . 


(B - 15) 


and the estimator Z has zero bias, i.e., it is unbiased. In the case of an unbiased 
estimator the covariance of the estimator is the same as the second moment of the 

A 

estimator about the true value. Let cov (Z) represent the covarance of Z. By 
multiplying both sides of Equation (B - 13) by B(Z) 
one obtains 

Z = Z + B(Z) A t (Z)Q- 1 (Y - X(Z ) (B - 16) 
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But Y = X<Z) and Y 


and 


By applying the "cov" 


cov (Z 


But 


cov (Z - 


Y = V, therefore 


Z-Z = B(Z)A x (Z)^ y ' 1 V , 


(B - 17) 


cov (Y - X(Z)) = cov (v) = Q y . (B - 18) 

operator to both sides of Equation (B - 17) we obtain 

Z)= l(Z)A T (Z)Q y _1 A(Z) |(Z) 

= B A t 2 " 1 AB 

= gr 1 g= S(Z). (B - 19) 

Z) = E [{(Z-Z) - E(Z-Z)}{(Z-Z) -E(Z-Z)} J 
= e[(z-E(Z))(z - E(Z)) t ] 

- cov (Z) 


64 



by definition where Equation (B - 15) was used and E(Z) = Z itself. Hence, 
Equation (B - 19) becomes 

cov (Z) = . B(Z) . (B - 20) 

According to Equation (B - 12) ^(Z) - £(Z). Using this approximation along 
with the definition of B, from Equation (B - 9), Equation (B - 20) becomes 

cov(Z) = [A T (Z) fi^ACZ)]" 1 = I(Z) • (B - 21) 

We wish to strongly emphasize that the usual assumption of the unbiasedness 
of the least squares estimate and the assumption that Equation (B - 21) represents 
the covariance of the least squares estimate rest on the linearity condition given 
by Equation (B - 12). The validity of this linearity condition is influenced by the 
degree of nonlinearity of the function which relates the state Z to the state of 
observations Y and by the distance between Z and Z. This last factor is highly 
correlated with v , the noise on the observations. 

Since ||(Z) is available at each step of the iteration procedure described by 
Equation (B - 11), then for the final ("best") estimate giving B(Z), assuming true 
convergence, Equation (B - 21) provides us with a means of calculating the 
covariance matrix of Z . This will be used in the next section (Appendix B.2) to 
obtain the error cone angle associated with the best estimate array Z. 
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For further detail concerning estimation theory in general see Deutsch 
(1965). 

B.2 The Construction of Error Cone Angles 

In Appendix B.l we provided a computational procedure for obtaining Z , 
the least squares estimate of a vector Z (and also X after the overdetermination 
equations are employed using the components of Z as the independent variables) . 
We also established the mean and variance of Z as a random variable by assuming 
unbiased measurements and a linearity condition. To generate error cones for 
n what we require as well is the precise distribution of the least squares estimate 
of B t and B 2 . To obtain this, one further assumption, previously stated but not 
used until now, is needed. That is, the noise v on the observations is assumed 
normally distributed. Notice that Equation (B - 17) gives the least squares 
estimate Z as 


, Z = Z + (a t (Z)fl y ” 1 4(Z))" 1 A T (Z)fi y - 1 v 


"V ^ 

where Z is, of course, a constant. Thus, Z is a linear function of v . A linear 
function of a normal random variable is normal also. A convenient feature of 
the normal distribution is that it is completely determined by its mean' and 
variance. The least squares estimator Z is unbiased. Hence, its expectation is 
Z. Its covariance matrix is given by Equation (B - 21). Thus, under the assump- 
tion of normal noise on the observations, the distribution of Z is completely 
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specified. However, our main interest now is focused on the distribution of the 

r\s r\j 

least squares estimates of the vectors Bj and B 2 , i. e. on the first six components 
of Z. 

Let Z r be the least squares estimate of and B 2 , i. e. Z = (!' , W y , W z ) 
[See Table 1 and Equation {m D - 5)] . Then, since Z is normal, Z' is also 

r\j r\j 

normally distributed with expectation the true values of B x and B 2 , The covariance 
matrix of Z can be obtained simply by deleting the last two rows and columns 
of the covariance matrix of Z (Graybill, 1961). Using this technique we can 
construct the exact distribution of the least squares estimate of B x and B 2 . Also 
with this background information we can perform useful simulation studies, as a 
test of the estimation scheme, and obtain least squares associated error cone 
estimates for simulation or real cases. 

The ultimate goal of our least squares error analysis is not to obtain best 
estimates of Bj and B 2 but to obtain the best possible estimate of the true shock 

r\j 

normal n. This is functionally related to B x andB 2 by 


rv, — 

n = F(B lt B 2 ) , 


(B - 22) 


where, from Equation (A - 6), 


n 


AB x (Bj xB 2 ) 

I AB x (B* xB 2 )| 


F(B 15 B 2 ) 


(B - 23) 
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with 


AB « B 2 ~ Bj . 

But since we do not possess the true values and B 2 (except in the eases of 
simulation studies), we must be satisfied with our least squares best estimates 
and B 2 and obtain an estimate n of n as 

n = F(B lf B 2 ) = F(Z') . (B - 24) 

Of major interest to us then is the statistical distribution of the error made in 
estimating n byn. A natural measure of this error is the angle between n and 
n . Thus, we define the function 0(n) as 

<£(n) a cos" 1 (n • ri) , (B - 25) 

where the principal value is understood. Obviously <£(n) is the angular error 
introduced by using n as an estimate of the true normal. Now we introduce 
another function \p( Z ' ) as 

0(Z') = 4>[F(Z')] . (B - 26) 
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0(Z' ) can be interpreted as the angular error in the normal estimate caused by 
using the least squares estimates of B t and B 2 rather than their inaccessible 
true values. In carrying out the least squares procedure for estimating normals, 
it is the statistical distribution of ’/'(Z' ) rather than that of Z' which is of in- 
terest. Since the function <j> is highly nonlinear the only reliable procedure for 
estimating the distribution of </<Z' ) is a Monte Carlo one. As constructed, the 
distribution of </»(Z' ) lies between zero and 180°, and our goal is to estimate a 
cone angle a. such that 95% of the distribution lies between zero and a. [The 
error cone geometrical interpretation is given in Section IV-A,] In other words, 
the probability is 0.95 that our estimate of the shock normal obtained by our least 

V 

squares procedure will lie in this cone. This 95% error cone clearly has in- 
tuitive appeal as a measure of our ability to estimate shock normals with the 
least squares procedure. It has the disadvantage, however, of being a single 
parameter measure of a cone that more precisely should not be described as 
being right circular. That is, strictly speaking the covariance matrix resulting 
from the least squares scheme contains enough information to be used to obtain 
a cone with an eliptical cross-section rather than a circular one. There is no 
reason to expect an "isotropic noise" situation to exist in general, and, in fact, 
there is good reason to expect otherwise for an average interplanetary type 
shock. But within the capability of the overall scheme, considering its limiting 
assumptions, the single parameter measure of an error cone should certainly 
be adequate. 
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With the covariance matrix of Z' in our possession we can implement the 
Monte Carlo process for estimating a by the following steps [Note: These same 
steps are used for obtaining the mean value error cone using its covariance 
matrix as discussed in Section IV-A] : 

1. Sample K times randomly from the distribution for Z ' as defined by its 
covariance matrix. 

2. Evaluate the function ^(Z' ) at each of these K points and thereby obtain 
K functional values of 0. 

3. Choose the functional value which is the smallest value that is larger 
than 95% of the K functional values obtained from step 2. This value is 
an unbiased estimate of a. 

The variance of the estimate of a obtained from steps 1, 2, and 3 above is 
inversely proportional to the Monte Carlo sample size K. From elementary 
probability theory it can be easily shown that, to a 95% certainty, the true 
percentage of the distribution of n contained in the 95% critical error cone 
centered at n is not -less than 94% if K > 2,420. For greater reliability we choose 
K = 3,000 for each cone angle calculation. 
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APPENDIX C 


COMPUTER PROGRAMS USED IN SCHEME 

The programs listed here, and used in this work, were written for 
the IBM 360-75 J computer. 

C.l Program to Obtain Best -Estimates ; Main and Subroutines 

On the eight pages following this preface there appears the shock program 
listing for real and simulated shocks. Table 5 shows the input format for the 
relevant input quantities. 


Table 5 
Input Format 



Quantity 

Format Designation 

Descriptive Notes 


'’iPRO 

1 1 

Integer with MFW* of 1 

Switches ^ 

ISWTCH 

I 1 

Integer with MFW of 1 


ICASE 

12 

Integer with MFW of 2, right adjusted 


^ ISAMPL 

I 2 

Integer with MFW of 2, right adjusted 

- 

XSTART 

8 F 6.2 

FPN** with MFW of 6 


SIG 

11 F 6.2 

FPN with MFW of 6 


NN 

1116 

Integer with MFW of 6, right adjusted 

- 

Y 

11 F 6.2 

FPN with MFW of 6 


' *MFW means maximum field width. 
**FPN means floating point number. 
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Recall that input XSTART, which is the starting vector for real shock studies, 
is Z t r ue (See Section IV-B) for simulated shock studies. [A sample of the printed 
output is given in Appendix D.l and described in Section IV-B.] The four switche 
comprise four numbers on a single input card as shown below: 


IPRO 


1SWTCH 


(CASE 


ISAMPL 


Switehl 
card J 


One digit 
per box 


The switch card is the first data card. The second data card is the XSTART 
card carrying eight numbers. The third card is SIG with eleven numbers. The 
fourth is NN with eleven numbers. And the last set of cards comprises the Y 
data array, each card of eleven numbers until all N ( = 2NN(i)J data points are 
listed. For example, if N is 28, then the first Y card has eleven numbers, the 
second eleven also, and the third card has six numbers. Table 5 shows this over- 
all order for the cards. The Y array is used only for real shock cases; for 
simulated shocks no Y data is necessary, of course. In all other respects the 
above comments hold for both real and simulated shock cases. If in the case of , 
a real shock a second or third, etc., XSTART is used, these are placed in order 
immediately after the Y array. And if another real shock is to be processed, 
the entire order of cards, from the 1st XSTART to the last XSTART, is repeated. 
And repeated again for a third shock, etc. But each' separate real shock case 
must have the same number of XSTART’ s. In the simulated cases of more than 
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one 'shock the set of XSTART, SIG, NN is simply repeated for each case. For 
any given computer run, for a real or simulated study, only one switch card is 
used and it is always the 1st card after the control cards. Table 6 schematically 
represents the card order for real and simulated cases. Below we describe the 
switches and how they are used. IPRO is used to control whether a real or 
simulated shock is to be pro cessed. ISWTCH can switch from XSTART unaltered 
to XSTART changed to XMEAN (from Y array) as the starting vector of the 


Table 6 

Input Data Card Order (One card for each line) 
Real Study Simulated Study 

Switch Card Switch Card 


XSTART (1) 

SIG 

NN' 

Y 


1st shock 
case 


XSTART 

SIG 

NN 




> 


J 



} 


1st shock 
case 


2nd shock 
case 


XSTART(2) 

XSTART(3) 

etc. 


( repeat order \ 
with same no. j 
of XSTART’ s / 


2nd shock 
case 
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iteration process. ICASE is just the number of different shocks (or cases) to be 


studied. ISAMPL is the number of samples of a given shock case (from the 
random number generator) to be considered for a simulated case, but it is the 
number of starting vectors for any given real case. [Note: Program is arranged 
to require the same number of "samples" for each shock case, for both real and 
simulated studies] . Table 7 shows what values the switches must have to perfori 
the duties described above. Section IV-B gives an example of their use for a 
simulation study. 


IPRO 


Table 7 

Program Switches 
( 0, real shock(s) being processed 
[_1, simulated shock(s) being processed- 


ISWTCH 


T 0, unaltered XSTART 

1_1, XSTART changed to XMEAN for iteration 


ICASE 

ISAMPL 


.from 1 to 99 equal to number of shocks 


from 1 to 99 equal to number of samples 
for each shock. 


74 



As an example of the use of the switches in a real shock study of two shocks 


using three starting vectors, the switches would be: 


IPRO = 0 1 


ISWTCH = 0 

ICASE = 2 


Switch card 


0 

0 

0 2 

0 3 


ISAM PL =3 J 


[Note: If ISAMPL is other than 1, ISWTCH should be 0 for a real case. If this 
is not adhered to, the program will successfully run but waste computer time 
through repeated operations] . 


For a real study of a single shock of N = 90 (number of components of Y ) 
with one starting vector the total program running time on the IBM 36 0-7 5 -J is 
only about 0.3 minutes. 
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uuuuuuuuuuuuuu 


IMPLICIT REAL48 (A-H,M,0-Z) 

DIMENSION X ( 15) ,G <7 ) , NN( 1 1 ) , Y ( 550 ) ,XXX( 8) ♦ XX(26,15), 
1B(8»8),A( 55 0*8) ,Q( 550 ),M(11),C(11,8),SIG(I1 ) , GK ( 1 1 ) 
DIMENSION XMEAN ( 15) »ZZ ( 16) , XSTART (15) , XL OSS < 1 6 ) * QUAL { 1 6 ) 
COMMON/XSQR/G4A,G5A,G6A,G4B, G5B » G6B , G4E » G5E , G6E * 

1 XN I NE » TEN, ELEVEN , XLOSSM , QUAL .ZZt QUALM » I S WTCH 
COMMON/AASUB/C 
CONV =57 • 2957 8D0 


IPRO TELLS THE PROGRAM WHETHER OR NOT WE WILL USE XMONTE TO 
GENERATE THE Y ARRAY (IPR0=1 USE XMONTE ( S I MULATED SHOCK) IPRO=0 
USE XMONTE (REAL SHOCK)) 


ISWTCH TELLS SUBROUTINE XSQ WHETHER' OR NOT TO SET XST ART =XME AN . I F 
ISWTCH=1 SET XSTART ( I ) ,FOR 1=1, 8, EQUAL TO XMEAN ( I ) IF ISWTCH=0 DO 
NOT SET THEM EQUAL. 

ICASE IS THE NUMBER OF DIFFERENT CASES TO BE PROCESSED 


IS AM PL FOR SIMULATED SHOCKS IS THE NUMBER OF SAMPLES PER CASE. FOR 
SHOCKS I SAM PL IS THE NUMBER OF XST ARTS PER CASE. 

RE AD (5, 09 ) I PRO, ISWTCH, I CASE, I SAMPL 

ISSAME=0 

ITOTAL=l 

I TOTLE — I CASE ^ I SAMPL 

8 FORMAT ( 1116) 

25 FORMAT ( 1H0 INPUT XSTART WAS * , 0 8 { F9 ♦ 3 , 1 X ) ) 

111 FORMAT ( 1H0, 'XSTART REPLACED BY XMEAN') 

113 FORMAT (1H , 'SIMULATED SHOCK BEING PROCESSED') 

112 FORMAT ( 1H , 'REAL SHOCK BEING PROCESSED') 

31 FORMAT ( 1H1 ,' SHOCK PROGRAM OUTPUT') 

6 FORMAT ( 11F6.2) 

27 FORMAT ( 1H0 ,' THE TOTAL NUMBER OF DATA POINTS, N, IS ',14) 

18 FORMAT ( 1H0, 'THE INPUT Y ARRAY WAS •) 

32 FORMAT ( 1H , 1 1 ( F9 . 3, IX ) ) 

9 FORMAT (21 1,212) 

S7 FORMAT ( 1H0 , 2 4X , 3HB1 X , 7X, 3HB1 Y , 7X , 3HB 1 Z , 7X . 3HB2X , 7X , 3HB 2Y , 

17X ,3HB2Z »7X ,3H WY,7X»3H WZ,7X»3H N1,7X,3H N2,7X,3H WX ) 

26 FORMAT ( 1H , ' THE INPUT SI G WAS * , 1 1 ( F9 . 3 , 1 X.) ) 

28 FORMAT (1H * ' THE INPUT NN WAS ',11 (19, IX)) 

55 FORMAT ( 1 HO ,' THE CORRESPONDI NG G VALUES ARE*) 

70 FORMAT ( 1H ,7<F9.3,1X)> 

56 FORMAT ( I H , 6X , 2HN1 , 8X , 2HN2 , 8X , 2HW X , 8X , 2HNX , 8X , 2HNY ,8X 
1 ,2HK'Z,8X,2H0P ) 

11 FORMAT ( 1H0 ,* THE NUMBER OF I TER AT I ONS , L , I S ',15) 

12 FORMAT ( 1H0, 'THE BEST ESTIMATE INDEPENDENT PARAMETER MATRIX IS') 

67 FORMAT ( 1H , 5X ♦ 3HB 1 X ,7X , 3HB 1 Y, 7X , 3HB1 Z , 7X ,3HB2X , 7X , 3HB2Y , 7X , 

1 3HB2Z »7X » 3H WY,7X,3H WZ,7X,3H N1,7X,3H N2»7X,3H WX ) 

50 FORMAT ( 1H , 84X , F 1 6. 5, ' =LOSS M'//) 

114 FORMAT (1H ,8 CF9 .3 , 1 X ) , 4X , FI 6 . 5 , ' =LOSS ',12,' Z= *,F8.4) 

44 FORM AT ( 1 H , 8 (F9 . 3 , 1 X ) , 4X , FI 6. 5 , *=LOSS ',12) 

106 FORMAT ( 1H0, 'THE BEST ESTIMATE DEPENDENT PARAMETER MATRIX IS ') 

1 07 FORM AT ( 1H , 16X, 2HN1 ,8X , 2HN2 , 8X , 2HWX , 8X , 2 HNX ,8X , 2HNY , 8X ,2HNZ. 
18X.2HOP) 

109 FORMAT ( 1H * 1 OX , 7 ( F9 . 3 , 1 X ) » 04X , F 1 6 . 5 , ' =QUAL ITY ',12) 

105 FORMAT ( 1H ,10HMEAN»S G'S, 7 ( F9. 3 , 1 X ) , 4X , F 1 6 .5 , ' =QUALITY M ' // ) 

16 FORMAT ( 1H0 , ' 8, THE COVARIANCE MATRIX OF FINAL ESTIMATE, IS •) 

29 FORMAT ( 1H , 1 1 (F9 .3 , IX ) , 'MEAN VALUES') 

40 FORM AT ( 1HO , * THE CONTRACTED FORM OF DERIVATIVE MATRIX, A, IS') 

23 FORMAT ( 1H , 8 (F9 .3, IX ) , 10X ) 

21 FORMAT ( 1H ,8(F9.5,1X)) 

71 FORMAT ( 1H0, 'AAVE=' , F7. 3-, 5X , ' ABE= • , F7 .3 , 5X, » AVE , AB E= • , F7. 3 ) 

1114 FORMAT (8F6. 2 ) 

1115 FORMAT ( 1H0, 'THE NUMBER OF DIFFERENT CASES TO BE PROCESSED IS ',13) 
3 READ (5,1114) (X(I ) *1=1*08) 

DO 115 1=1,8 

1 15 XXX< I )=X (I ) 

READ ( 5 , 6) SI G 
READ (5,8 )NN 

NTOT =NN ( 1 ) + NN ( 2 ) +NN ( 3 ) +NN ( 4 ) +NN ( 5 ) + NN (6) +NN ( 7 ) +NN ( 8 ) +NN ( 9 ) +NN ( 10) 

1 +NN ( 1 1 ) 

230 WRITE (6,31) 

1 F ( I SWTCH • EQ • 1 ) WR I TE ( 6 » 1 1 1 ) 

IF(IPRO.EQ.O)WRITE(6,112) ' 

I F ( I PRO «EQ . 1 ) WR ITE(6» 113) 

WRITE (6, 111 5) ICASE 

WRITE (6, 27) NTOT 

IF ( IPRO.EQ. 1 )G0 TO 24 
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IF ( ( IPRO.EQ.O ) . AND. ( I SSAME. NE.O ) ) GO TO 120 
DO 17 1=1 iNTOT ill 
IKK = I +10 

17 READ(5,6) (Y(JXC>, JXC=I ,IKK) 

120 WRITE ( 6 » 18 ) 

DO 19 I =1 jNTQT , 1 1 
I KK = 1 + 10 

WRITE (6, 321 (Y ( JXC ) , JXC=I , I KK ) 

19 CONTINUE 

24 WRITE (6, 57) 

WRITE (6) 26 ) S I G 
WRITE (6 (28 ) NN 
CALL CON (X )G ) 

WRITE (6, 25) (X (I ) ,1=1,08) 

WRITE (6,55) 

WRITE (6 ,56) 

WRITE (6,70 )G 

IF ( IPRO »EQ » 1 ) CALL XMONTE ( SI G , NN ♦ Y , X ) 

CALL XSQ (SIG ,NN,Y ,X , XX , L , B » XLOSS , XMEAN ) 

WRITE (6 , 1 1 ) L 
V'RITE (6,12) 

WRITE ( 6 ,67) 

JT0P=L+1 

DO 13 1=1, JTOP 

IF (I .EO. 1 ) WRITE (6 ,29) ( XMEANCLM ) ,LM=1 ,8) , XNI NE , TEN , ELEVEN 
IFd.EQ.l )WRITE (6, 50 )XLOSSM 
NC OU NT = 1 — 1 
J JTOP= 1 6 

IF (I .EO. J JTOP) WRITE (6,44) ( XX (I , LM ) ,LM=1 ,8) , XLOSS (I ) ,NCOUNT 
IF ( I .NE . JjTOP ) V/RITE (6, 11 A) (XX (I »LM) , LM = 1 ,8) , XLOSS ( I ) ,NCOUNT,ZZ ( I ) 
IF (I .EQ . JTOP )G0 TO 100 
GO TO 13 

100 DO 101 11=1,8 

GK (II )=XX (I » I I ) 

101 CONTINUE 
13 CONTINUE 

WRITE (6,106) 

WRITE (6, 107) 

DO 102 1=1. JTOP 
NC OUNT = 1 — 1 

IF (I .EQ. 1 ) WRITE <6 , 105) (XMEAN(LM) ,LM=9» 15) , QUALM 
WRITE (6, 109 ) (XX ( I , LM ) ,LM=9 , 15 ) , QUAL ( I ) .NCOUNT 

102 CONTINUE 

CALL AA (GK , NN , A ) 

WRITE (6, 16 ) 

DO 41 1=1,8 

WRITE (6.21 ) ( BC I ,LM) ,LM=1 ,8) 

41 CONTINUE 

WRITE (6,40 ) 

DO 20 1=1,11 

WRITE (6,23 ) ( C(I ,LM) , LM=l ,8) 

20 CONTINUE 

CEA=( (G4E*G4A)+(G5E*G5A)+(G6E*G6A) ) 

CEB=( (G4E*G4B) + (G5E*G5B);MG6E*G6B) ) 

CAB=( (G4A*G4B)+(G5A*G5B)+CG6A*G66) ) 

CC=DARCOS(CEA) 

BB=DARC0S( CEB) 

DD=DARC0S ( CAB ) 

AAVE=CCSCONV 
ABE=BB#CONV 
AVEABE=DD*CONV 
WR I TE (6.71 ) AAVE ♦ ABE, AVEABE 
IF ( I T OTAL. EQ • I TOTLE )GO TO 99 
ITOTAL=ITOTAL+l 
ISSAME=I SSAME+1 
IF (I SSAME .EQ . I SAMPL ) I SSAME=0 
I F ( IPRO • EQ .0 ) GO TO 118 
DO 116 1=1,8 
116 X(I )=XXX (I > 

GO TO 119 

118 I F ( (IPRO.EQ.O).AND. ( I SSAME. NE . 0 ) ) READ < 5 , 1 1 1 4 ) ( X ( I ) , I =1 , 08 ) 

119 IF ( ISSAME .EQ.O )G0 TO 3 
GO TO 230 

99 STOP 

END 

SUBROUTINE C0N(X,G) 

IMPLICIT REAL£*8 (A-H.M.O-Z) 

REAL*8 Ml . N 2 
DIMENSION X(15)»G(7) 
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bix=x< n 

B 1Y=X (2) 

B 1 Z =X ( 3 ) 

B2X=X (4) 

B2Y =X ( 5 ) 

B2Z=X (6) 

WY=X(7) 

WZ=X (8 ) 

XX=B2X— BIX 
Y=S2Y-B1 Y 
Z=82Z— B1 Z 
PIE=3. 1415927 

F=(XX**2+Y**2+Z**2J/(4.0D0*PIE> 

A=B2X**2+B2Y S42+B2Z **2-8 1 X*#2-B1 Y*£2-B1 Z**2 
QX= (BlY*B2Z )-(BlZ*B2Y ) 

QY=(B1Z*B2X)-(B1X*B2Z ) 

0Z=(B1X*B2Y )-(BlY*B2X ) 

MX=(QZSY )— (QY^Z ) 

MY = (QX*Z )-(QZ*XX ) 

MZ=( QY*XX ) - ( QX*Y ) 

E=- ( <WY*QY ) +( WZ*QZ ) ) XQX 
T=E*XX+WY*Y+WZ*Z 
D=E*MX+WY *MY+WZ*MZ 
6G=B1 X*MX+B1Y*MY+B12*MZ 
M = DSQRT (MX*=S=2 + MY**2-S-MZ**2 ) 

SX=WY*QZ-WZ*QY 
SY=VJZ*QX-E*QZ 
SZ=E *QY— \VY Wx 

R= ( 62X*SX+B2Y 4SY+B2Z4SZ ) / C B1 X4SX+B1 Y4SY+B1 Z*SZ ) 

RR=— ( (F4GG ) / ( T£D ) ) 

WX=E 

Ml= ( (R-l .000 ) /R ) 4 ( 597 9.1 4D0#RR ) 

N2= (R-l .ODO )£ ( 5979 .1400 tRR ) 

G( 1 )=N1 
GC2 )=N2 
G(3)=WX 
G ( 4 ) -MX /M 
G (5 ) =HY /M 
G ( 6 ) ~MZ /M 

G (7 ) = ( ( ( RR — D**2 ) / ( M4. 2 ) ) -( A/ ( 8 • *P I E ) ) ) 

RETURN 

END 

SUBROUTINE FI(X.NN.F) 

IMPLICIT REALMS ( A— H . 0— Z ) 

DIMENSION X ( 1 5 ) » NN (11) , F ( 550 ).M( 11 ) »G(7) 

C FI CALCULATES EXPECTED VALUE OF OBSERVATIONS AND STORES THEM IN F 

C NN(I) IS NUMBER OF OBSERVATIONS OF ITH VARIABLE 

C X ARRAY CONTAINS SHOCK PARAMETERS 

CALL CON(X.G) 

X(9)=g(1 ) 

X( 10)=G(2> 

X ( 1 1 ) =G ( 3 ) 

N = 0 

DO 1 J=1 , 1 1 

N=N+NN( J ) 

1 M(J)=N 
J=1 

DO 2 1=1 ,N 

IF (I.LE.M(J))GO TO 2 
J= J+ 1 

2 F ( I ) =X { J ) 

RETURN 

END ’ 

FUNCTION XL(X.SIG.Y.NN) 

IMPLICIT REALMS (A-H.O-Z) 

DIMENSION X(15)»SIG(11),Y( 550 ) , Q( 550) ,NN(11),B(11) ,F (550 ) ,M ( 1 1) 

1 .G (7) 

C XL IS THE LOSS FUNCTION 

C X IS THE ESTIMATE OF SHOCK PARAMETERS 

C SIGtl) IS THE SIGMA VALUE ON ITH VARIABLE 

C Y(I> IS THE ITH OBSERVATION 

C NN ( I ) IS NUMBER OF OBSERVATIONS ON ITH VARIABLE 

J=1 
N=0 

DO 1 1=1*11 

8(1 )=1.0D0/(SIG(I >*SIG(I )) 

N=N+NN(I ) 

1 M(I)=N 

DO 2 1=1. N 
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IF(I .LE.M(J) ) GO TO 2 
J=J+1 

2 QCI )=0( J) 

CALL FI (X,NN,F) 

DO 3 I =1 ,N 

3 F { I )=F { I )-Y (I ) 

XL=0 .ODO 

DO 4 I = 1 » N 

4 XL=XL+(Q (I )*FU )*F(I ) ) 

RETURN 

END 

SUBROUTINE AA ( X » NN , A ) 

IMPLICIT REALMS (A-H,0-Z) 

DIMENSION X ( 1 5) , NN{ 1 1 ) , A < 550 ,8),B(6),ID(8,8),XK(3,6) 

1 <C( 1 1 >8) »YM(3,6)»D1 (6)»G1(6),T1(6),F1 (6) »XJ(6)»SS(3,6) » XET A ( 6 ) » 

2 XX I (6) ,W(3) ,G (3) ,XM(3) , S(3) ,H(6> , WS< 3,3) ,WE(3 ) , WW ( 3 ) ,V<3) . TVJ(3) , 
3DW(3 ),U(3)»M(11> 

COMMON/A ASUB/C 
DO 1 1=1*6 

1 B < I ) =X (I ) 

W(2)=X (7 ) 

W ( 3 )=X (8) 

XX=B{4 )— B( 1 ) 

Y=B (5)-B( 2) 

Z=B ( 6 ) — B ( 3 ) 

PIE=3. 14159 27D0 

F=(XX**2+Y**2+Z**2) /(4.0DO*PIE) 

Q(1)=(B(2)*B(6))-(B(3)*B(5>) 

Q(2)={8(3)*B(4))-<B(1>*BC6>) 

Q(3) = (B(l)tB(5))-(B(2) *B(*4) ) 

XM( 1 )=(Q (3)*Y)-<Q<2)*Z ) 

XM(2)= (Q( 1 ) *Z )-(Q(3)*XX) 

XM(3 )= (Q (2 )*XX )-(Q ( 1 >*Y) 

E=-( <W(2)*Q (2))+(W(3)*Ol3) ) )/Q< 1 ) 

T=E*XX+W (2 )*Y+W( 3)*Z 

D=E*XM ( 1 )+W (2)*XM(2)+W(3)*XM(3) 

6=B ( 1 ) *XM ( 1 )+B ( 2)*XM(2 )+B(3 )*XM( 3) 

S ( 1 ) =W ( 2 ) *0 ( 3 ) -V/ ( 3 ) SQ ( 2 ) 

S(2)=W(3)*Q{ 1)-E*Q<3> 

S ( 3 ) =E *Q ( 2 ) — W ( 2 ) *0 ( 1 ) 

ETA=B(4)*S( 1 )+B(5)=rS{2)+B(6)*S(3) 

X I =B( 1 )*S(1 )+B(2)*S(2)+B(3)*S(3) 

R=ET A/XI 

RR=-( (F*G >/<T*D) ) 

DEN= (R— 1 .ODO )*(5979.14D0#RR ) 

DO 2 1=1*8 
DO 2 J=1 >8 
10 (I *J )=0 

I F ( I .EQ . J ) I D ( I * J) = l 

2 CONTINUE 
DO 3 1=1*6 

XK ( 1*1 )=B (6 )*I DC2, I )+B(2)*ID(6, I )-B( 5)*ID(3,I ) 

1 -B(3)*ID(5,I) 

XKC2*I )=8(3)*ID(4,n+B(4)*ID(3,I)-B( 1)* ID (6,1 ) 

1-6(6)4=10(1.1 ) 

XK (3,1 )=B(1)*IDC5,I)+B(5)=5=ID(1,I)-B(2)*ID(4,I) 

1-B (4)*ID(2,I ) 

3 CONTINUE 
DO 4 1=1,6 

C ( 1 1 »I )=- ( ( E=v=XK ( 1 , I )+W( 2>*XK<2, I)+W(3>=4=XK(3»I))/ r (Q(l))) 

4 CONTINUE 

. DO 5 1=2,3 

5 C(Il,I+5)=— ( (Q(2)*ID(2,I )+G(3)*ID(3,I))/(Q(l)) ) 

DO 6 1=1,6 

YMC 1 , I )=XK(3 ,1 ) 58= Y — XK (2,1 )*Z+Q( 3>*( ID<5,I)-ID(2,I)) 

1-QC2)*(ID(6,I )-ID(3,I ) ) 

YM ( 2 , I ) =XK ( 1 , 1 )*Z-XK(3,I )*XX+Q( 1 )*( ID (6 , I ) - ID ( 3 , I) ) 

1-Q (3)* (ID (4,1 )-ID( 1 , I ) ) 

YM (3,1 )=XK(2 , I )*XX-XK( 1 , I )*Y+Q (2)*( ID ( 4 , I ) - ID ( 1 , I ) > 

1—0 ( 1 ) * ( I D (5,1 )— ID(2»I ) ) 

OKI )=E*YM (1 , 1 )+W (2)*YM(2, I )+W(3 )*YM(3, I )+XM ( 1 )*C (11,1 ) 

G 1 ( I )=O.ODO 
DO 7 J= 1,3 

7 G1 ( I )=G1 ( I )+ID( J, I )*XM( J)+B{ J)*YM( J, I ) 

T1(I )=E4=(1D(4, I )-ID( 1 , I ) )+W(2) *C ID( 5, I ) — ID (2, I ) >+W(3 ) 
1*(ID(6,I)-ID(3,I) )+XX*C( 1 1,1) 

F1(I )={XX*(ID(4,I ) — ID( 1,1 ) )+Y* (ID(5,I)-ID(2,I))+Z 
1#(ID(6,I)-ID(3,I)))/(2.0D04=PIE) 

X Jd )=RR*(F1 (I ) /F+Gl ( I )/G-Tl (I ) /T-Dl ( I ) /D ) 
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SS { 1 , 1 )=W (2 )*XK< 3, I )-W (3)4XK<2, I ) 

SS (2,1 )=W(3> *XK( 1 , I )-E*XK<3, I > -Q < 3 ) *C < 1 1 » I > 

SS (3,1 ) =E=S=XK ( 2,1 )-W( 2)4XK (1,I)+Q(2)*C(11»I> 

XET A ( I )=0.0D0 
XXI (I )=0 .ODO 
DO 8 J=1 ,3 

XXI (I )=XXI ( I >+ID( J, I )*S( J )+B( J)*SS( J, I ) 

8 XET A (I )=XETA ( I )+ID( J+3, I )4S ( J) +B( J+3 >*SS( J > I ) 

H ( I )=(XI*XETA(I )-ETA*XXI (I ))/<XI**2) 

C(10»I )=DEN4(H(I )/(R-l .ODO)+XJ ( I )/RR) 

C ( 9 » I >=(R*C(10»I )-DEN*H(I ))/<R4*2) 

6 CONTINUE 
DO 9 1=2,3 

WS(1,I)=Q<3)*ID(2,I)-Q(2)#ID(3,I) 

WS (2, I )=Q (1 )*ID( 3, I )— Q (3)*C( 11 » 1+5) 

WS 13 » I )=Q ( 2 ) ( 1 1 » I + 5)— Q ( 1 ) *1 D ( 2 , I ) 

WW(I ) =0 . ODO 
WE ( I ) =0 ♦ ODO 
DO 10 J= 1,3 

WE (I )=WE (I )+B( J+3)*WS( J, I ) 

1 0 WW ( I ) =WW { I )+B ( J ) *WS ( J » I ) 

V(I ) = { X I #W E ( I )-ETA*WW(I ) )/(XI**2) 

TVMI )=XX*C< 11 ,I+5)+Y*I0(2,I )+Z*ID(3, I ) 

DW(I )=XM (I )*C( ll,I+5)+XM(2)*ID(2, I ) +XM < 3 ) * ID (3 , 1 ) 

U(I )=-RR*(TW(I ) /T+DW ( I )/D) 

C(10,I+5)=DEN*(U(I ) /RR + V ( I )/<R-1.0DO) ) 

C ('9 ,1+5 >= CR*C{ 10 ,I+5)-DEN4V( I ) )/(R*42) 

-9 CONTINUE 
DO 11 1=1,8 

DO 12 J= 1 * 8 
CC J,I )=ID(J,I ) 

12 CONTINUE 

11 CONTINUE 
N=0 

DO 13 L= 1 » 1 1 
N=N+NN (L ) 

13 M(L)=N 

DO 15 K= 1,8 
J=1 

DO 14 I=1,N 

IF < I .LE.M ( J ) >GO TO 14 

J=J+ 1 

14 A ( I ,K)=C( J,K> 

15 CONTINUE 
RETURN 

END __ 

SUBROUTINE XSQ (SIG » NN , Y , XSTART , XX , L, B, XL OSS , XMEAN ) 

IMPLICIT REALMS (A-H,0-2) 

DIMENSION SIG( 11 ), NN (11 ) ,Y(550 ) .XSTART ( IS) , XX ( 26 , 15 ) , ZZ ( 16 ) , 
18(8,8 ) ,XLOSS( 16) , XMEAN ( 15) , Q ( 550 ),G(7),G1(7),C(8),M(11) ,QUAL( 16) 
COMMON /XSG)R/G4A,G5 A ,G6A»G4B» G5B, G6B, G4E,G5E» G6E , 

1XNINE , TEN, ELEVEN , XLOSSM , QUAL , 2 Z , QUALM , ISWTCH 
SIG(I) IS SIGMA VALUE ON ITH PARAMETER 
NN(I) IS NUMBER OF OBSERVATIONS OF ITH PARAMETER 
Y ( I ) IS ITH OBSERVATION, OR ITH COMPONENT OF DATA ARRAY 
XSTART IS STARTING VECTOR FOR NUMERICAL SOLUTION 
XX(I,J) IS ITH ESTIMATE OF JTH PARAMETER 
L IS NUMBER OF ITERATIONS USED 

B IS FINAL PROPAGATED COVARIANCE MATRIX OF BEST ESTIMATE 
XLOSS(I) IS VALUE OF LOSS FUNCTION FOR ITH STEP OF ITERATION 
XMEAN IS DATA AVERAGE VALUE PARAMETER ARRAY 
MTT =0 
N=0 

CNORM =0 » ODO 
J = 1 

XN0RM=0 .ODO 
DO 1 1=1,11 
XMEAN (I )=0 .ODO 
N =N+NN (I ) 

1 M(I)=N 

DO 2 1=1, N 

IF ( I *LE • M ( J ) ) GO fO 3 

J=J+1 

3 Q (I ) = 1.0D0/< SIG ( J )*SIG ( J) ) 

2 XMEAN ( J ) =XMEAN (J)+(Y(I )/NN(J)) 

XNINE=XMEAN(9) 

TEN=XMEAN(10) 

ELEVEN=XMEAN (11) 

XLO SSM=XL (XMEAN, SI G ,Y,NN) 
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QUALM =DSQRT< N/XLOSSM) 

CALL CON C XMEAN * G ) 

G4A=G (4) 

G5A=G (5) 

G6A=G (6 ) 

CALL C ON (XSTART $G 1 ) 

G4E=G 1(4) 

G5E=G 1(5) 

G6E=G 1(6) 

DO 4 1=1,7 
J = I +8 

XMEAN ( J)=G(I ) 

4 XST ART ( J )=G1 (I ) 

IF ( 1SWTCH .EG. 1 )GO TO 10 
GO TO 11 

10 DO 12 1=1, IS 

12 XST ART ( I )=XMEAN < I ) 

11 DO 5 1=1,15 

5 XX (1,1 )=XSTART( I ) 

XLOSS ( 1 ) =XL (XSTART , SIG , Y » NN ) 

QUAL ( 1 ) = DSQR T ( N/X LOSS ( 1 ) ) 

CALL P ( XSTART , 0 , NN , V , B , C ) 

DO 20 L = 2, 16 
DO 6 1=1,8 

CNORM=CNORM+ (C ( I ) =XC < I ) ) 

6 XNORM =X NO RM+ ( XSTAR T ( I ) * XSTART ( I ) ) 

XNORM=DSQRT( XNORM ) 

CNORM =DSQRT ( CNORM ) 

2=CN0'RM/XN0RM 

IF ( L .EG • 2 )ZZ ( 1 )=Z 

IF(Z .LE. .0 10D0 )G0 TO 25 

mtt=l~i 

ZZ (MTT >=Z 
DO 7 1=1,8 

7 XSTART (I )=XSTART ( I ) -C ( I ) 

XLOSS ( L )=XL( XSTART, SIG ,Y f NN) 

QUAL ( L ) =DSQR T (N/X LOSS ( L ) ) 

CALL CON (XSTART ,G ) 

DO 8 1=1 ,7 
J = I+8 

8 XSTART (J )=G< I ) 

DO 9 1=1 , 15 

9 XX (L, X )=XSTART(I ) 

CALL P (XSTART ,Q , NN * Y ,8,0 
20 CONTINUE 
25 KMT =MTT+1 
ZZ (KMT )=Z 
L = MTT 

CALL CON (XST ART , G ) 

G4B=G (4 ) 

G5B=G (5) 

G6B =G (6 ) 

RETURN 

END 

SUBROUTINE P (X , G , NN , Y , B , C ) 

IMPLICIT REAL*8 ( A-H , 0-Z ) 

REAL #4 EPS 

DIMENSION U ( 8 * 8 ) , AFLAG (8) ,ATEMP(8) , Q ( 550 ) . X ( 1 5 ) , NN ( 1 1 ) , Y (550 ) , 
1 0(8,550 ) ,F(550 ) , A < 550 , 8 ) ,B(8,8),C (8) ,G(7) 

NC = 8 
MR =8 
NR = 8 
EPS=3 .O 
N = 0 

CALL AA ( X , NN , A ) 

DO 1 1=1,8 
C (I )=0 .000 
DO 2 J=l,8 

2 B(I,J)=O.ODO 
1 CONTINUE 

DO 3 1=1,11 

3 N=N+NN (I ) 

DO 4 1=1 ,8 
DO 5 J=1 ,8 
DO 6 L= 1 , N 

6 B(I ,J)=B(I , J)+(A(L,I )-0(L)*A(L, J) ) 

5 CONTINUE 

4 CONTINUE 

CALL GINV2 (B,U, AFLAG , ATEMP , MR , MR , NC * NR A NK , EPS ) 
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DO 7 1=1,8 
DO 8 J=1 ,8 

8 B ( I i J) =B ( J , I ) 

7 CONTINUE 

CALL FI (X »NN»F ) 

DO 9 1=1 ,N 

9 F ( I ) =F ( I ) -Y < I > 

DO 10 1=1,8 

DO 11 J=1,N 

11 D ( I , J )=0 »0D0 
10 CONTINUE 

DO 12 1=1,8 
DO 13 J=1»N 
DO 14 L= 1,8 

14 D{ I , J )=D (I , J) + (B(I ,L )*A ( J,L) ) 

13 CONTINUE 

12 CONTINUE 

DO 15 1=1,8 
DO 16 J = 1 » N 

16 C(I)=C(I)+(D(I,J)*Q<J)*F(J)) 

15 CONTINUE 
RETURN 

END 

SUBROUTINE X MONTE ( SI G , NN , Y , X ) 

IMPLICIT REAL *8 ( A-H , 0-2 ) 

DIMENSI ON SIGUl ),NN(11 ) , Y { 550 ),X(15),M(11),G(7> 
C SIG ARE SIGMA VALUES 

C NN GIVES NUMBER OF EACH TYPE OF MEASUREMENT 

C X IS TRUE SHOCK PARAMETERS 

C Y IS MONTE CARLO SAMPLE OF MEASUREMENT 

N=0 
J = 1 

DO 1 1=1,11 
n=n+nn ( I ) 

1 M ( I >=N 

CALL C ON { X , G ) 

X ( 9 ) =G ( 11 

X(10 )=G(2> 

XI 11)=G(3) 

DO 2 I = 1»N 

IF (I ,LE.M(J ) ) GO TO 3 
J=J+1 

3 Y ( I )=X { Jl+BARNl ( -1 , 1,1 2787, S I G{ J ) ) 

2 CONTINUE 
RETURN 

END _ , 

SUBROUTINE GAUSS ( KI X /, S , AM , V , H ) 

IMPLICIT REALM'S (A-H, 0-2) 

K “H 

A=0 • ODO 

DO 50 1=1, K 

CALL RAMDUdX ,IY ,V) 

I X = I Y 
50 A=A4-Y 

H0=H/12. 

H2=H/2. 

V= ( S* ( A-H2 ) ) /DSQRT ( HD ) +AM 
RETURN 

END 

SUBROUTINE R ANDU ( SI X/ , I Y , Y FL ) 

IMPLICIT REAL*8 (A-hUO-Z) 

DATA JJJ5/1027/ 

I Y=IX*JJJ5 
IF<IY)5,6,6 

5 I Y = I Y + 2 147483647 + 1 

6 YFL= I Y 

YFL=YFL*. 46566 13D-9 
RETURN 

BIO 

FUNCT ION BARN1 ( I , IKEY , IFRN.SD) 

IMPLICIT REALMS (A~H*G~Z) 


C SD THE DESIRED STANDARD DEVIATION 

C AMEAN THE DESIRED MEAN 

C H THE POPULATION SIZE 


DATA AMEAN/O • DO / 
DATA I HERE / 127 87 / 
DATA H/36.D0/ 

I F { I KEY ) 5 * 4 1 4 
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4 I HERE = I FRN 

5 IF ( I >6 *7 .7 

6 CALL GAUSS ( I HERE, SD , AMEA N, VAL,H) 

IFRN=IHERE 

GO TO 8 

7 CALL RANDUUHERE, IFRN, VAL ) 

I HE RE = I FRN 

8 BARN 1 =VAL 
RETURN 

END 

SUBROUTINE GINV2(A,U, A FLAG , ATEMP * MR * NR * NC , NR1 ♦EPS) 
DOUBLE PRECISION FAC * DOT , DOT 1 * DOT2 , TOL ,DSQRT 
DOUBLE PRECISION A ( MR , NC ) » U ( NC » NC ) * AFLAG ( NC ) , ATEMP (NC ) 
DO 10 I =1-9 NC 
DO 5 J=1 9 NC 

5 U ( I , J)=0 . 

10 U (I * I ) = 1 ♦ 

FAC =DOT ( MR ♦ NR , A » 1* 1 ) 

F AC = 1 • /OSQRT CFAC) 

DO 15 1=1, NR 
15 A(I , 1 )=A Cl , 1 ) ❖FAC 
DO 20 1=1, NC 
20 U(I , 1 )=U(I , 1 ) £FAC 
AFLAG ( 1 ) = 1. 

N=56 
NR 1 =NC 

TOL=( 10***EPS*,5**N)**2 
DO 100 J=2,NC 
D0T1=D0T(MR,NR , A', J, J) 

JM 1 = J— 1 
DO 50 L= 1,2 
DO 30 K=1,JM1 

3 0 ATEMP ( K) =DOT(MR , NR, A, J,K) 

DO 45 K=1 , JM1 
DO 35 1=1, NR 

35 A ( I , J)=A( I » J ) -ATEMP { K ) * A ( I ,K>*AFLAG<K) 

DO 40 1=1, NC 

40 U ( I , J) =U( I , J)-ATEMP <K)*U(I ,K> 

45 CONTINUE 
50 CONTINUE 

D0T2=D0T ( MR , NR ♦ A , J , J ) 

IF ( (DOT2/DOT1 ) -TOL) 55, 55, 70 
55 DO 60 1=1, JM1 
ATEMP (I >=0 * 

DO 60 K= 1,1 

60 AT EMP ( I ) = ATEMP (I )+U(K, I )*U(K, J ) 

DO 65 I =1 , NR 
A ( I , J ) =0 ♦ 

DO 65 K=1 , JM1 

65 A (I , J)=A( I , J)-A( I »K)£ATEMP(K }# AFLAG ( K) 

AFLAG ( J)=0. 

FAC =DOT CNC , NC , U, U , J ) 

F AC = 1 ,/D SORT (FAC) 

NR 1 =NR 1-1 
GO TO 75 
70 AFL AG ( J ) = 1 • 

F AC= 1 ,/DSORT (D0T2 ) 

75 DO 80 1=1 ,NR 
80 A(I , J)=A(I , J)*FAC 
DO 85 1=1, NC 
85 U(I , J)=U(I , J )*FAC 
100 CONTINUE 

DO 130 J= 1 , NC 
DO 130 1=1, NR 
FAC = 0 . 

DO 120 K= J 9 NC 
120 FAC=FAC+A( I , K)*U( J, K) 

130 A ( I , J)=FAC 
RETURN 

END 

FUNCTION DOT ( MR ♦ NR ♦ A, J,K) 

DOUBLE PRECISION A ( MR , 1 ) , X , DOT 
X =0 • DO 

DO 50 1=1 ,NR 
X =X + A ( I , J )*A(I ,K) 

50 CONTINUE 
DOT =X 
RETURN 
END 
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C.2 Program to Generate Error Cone Angles; Main and Subroutines 


On the five pages following this preface there appears the error cone program 
listing for real and simulated shocks. Table 8 shows the input format for the 
relevant input quantities. X refers to the first six components of Z true for a 
simulation study or of either Z mean or Z (best estimate) for a real shock study 
(program must be run twice for real cases — See Section IV- A). SIG and NN 
refer to the first six components of. the SIG's and NN's corresponding to the 
associated shock program, and N is the Monte Carlo sample size number (K in 
Appendix B.2) which is usually set equal to 3,000. B1 is the matrix B(Z) with 
the last two rows and columns deleted\ (See Appendix B.2 for explanation). The 
first four rows of Table 8 represent the first four data cards of the program in 
the order shown. The next six cards are the next six rows of the B1 matrix, 
respectively. 


Table 8 

Input Data Format 


Quantity 

Format Designation 

Descriptive Notes 

N 

I 4 

Integer with MFW* of 4, right adjusted 

X 

6 F 5.2 

FPN** with MFW of 5 

NN 

6 I 3 

Integer with MFW of 3, right adjusted 

SIG 

6 F 5.2 

/every \ 

FPN with MFW of 5 . 

B1 

6 F 9.5 1 row of 1 
\ matrix j 

FPN with MFW of 9 


*MFW means maximum field width. 

* *FPN means floating point number. 
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The printed output of the program is given, in sample form, in Appendix 
D.2 and is described in Section IV-B. The total program running time for 
N = 3,000 and 


L 

5=1 


NN(i) 


90 


on the IBM 360-75 J computer is only about 3 minutes. 



nonnnooo 


IMPLICIT REAL £8 (A-H.O-Z) 

DIMENSION B1(6,6),SIG(6),NN(6),X(6) 

READ ( 5 »8 ) N 
READ (5*1 )X 
READ(5,2)NN 
READ (5.1 )SIG 
WR I TE ( 6 » 3 ) X 
WRITE (6 » 6 ) SI G 
WRITE (6,4 ) N 
WRITE (6.5) NN 

5 FORMAT (1H , • THE INPUT NN VALUES WERE ',6(15* IX)) 

3 FORM AT (1H1, 'THE INPUT X VALUES WERE • ,6 ( F5.2 . 1 X ) ) 

4 FORMAT (1H .'THE INPUT VALUE FOR N WAS «,I5) • 

6 FORMAT ( 1 H , * THE INPUT SIG VALUES WERE • , 6 ( F5 .2 , 1 X ) ) 

WR I TE (6,7 ) 

7 FORMAT ( 1HO THE INPUT VALUES FOR B1 WERE •) 

DO I 6 1=1.6 

11 FORMAT (6F9. 5 ) 

READ (5. 11 ) (B1 (I , J) , J=1 ,6) 

10 WRITE (6. 12) ( B1 ( I , J) , J=1 ,6) 

12 FORMAT ( 1H . 6(F9«5,1X)) 

1 FORMAT <6F5. 2) 

2 FORMAT <6I3) 

8 FORMAT (14) 

CALL CONE (Bl, SIG, NN.X.N) 

RETURN 

END 

SUBROUTINE F IND ( MM » A ) 

IMPLICIT REAL*8 (A-A.O-Z) 

DIMENSION A(IOOOO) 

XM=MM 

ISIG1=XM*( .05D0 ) + 1 . 0D0 
IS IG2=XM* ( . 0 1 DO )+l *000 
ISIG3=XM4( .005D0 )+l .ODO 
WR I TE ( 6 ,900 )ISIG1,ISIG2»ISIG3 
900 FORMAT ( 1H0 , • ISIG1= ',13, » ISIG2= ',13,' ISIG3= ',13) 

DO 15 I=1,ISIG1 
X=AC 1 ) 

K=1 

DO 5 J=2,MM 

IF (X.GT.A(J) )G0 TO 5 

K=J 

X = A( J) 

5 CONTINUE 
A(K)=-10.0DO 

IF (I *EQ*ISIG1.0R*I • EQ » I SI G2 • OR ,I,EQ»ISIG3)WRITE(6, 1969) I, X 
1969 FORMAT ( 1H0, 'THE *,I3,' VALUE VMS '.D17.8) 

15 CONTINUE 
RETURN 

END 

SUBROUTINE CONE ( Bl , S IG , NN , X , N ) 

IMPLICIT REALS8 (A-H.O-Z) 

DIMENSION Bl{6,6)'» SIG(6),NN(6),X(6),XNN(6) ,A1 (6),Z1(6) 

1, T2(6),E1(6),T1 (6)»FI1 (10000) , FI 2 ( 1 0000 ) »XM 1 (6 ) ,XM2(6) 

2 , XM (6 ) , B2 (6 ) 

Bl IS COVARIANCE OF L.S. ESTIMATE 
SIG IS DEVIATION ARRAY 

NN ARRAY GIVES NUMBER OF READINGS IN EACH DATA CHANNEL 
X ARRAY GIVES TRUE SHOCK PARAMETERS 
N IS MONTECARLO SAMPLE SIZE 
SUBSCRIPT 1 REFERS TO L.S. ESTIMATE 
SUBSCRIPT 2 REFERS TO MEAN VALUE ESTIMATE 
DATA ON L.S. CONE IS PRINTED FIRST 
CALL XNORM (X »XM ) 

DO 1 1=1,6 
XNN Cl )=NN (I ) 

1 B2(I ) = SI G (I ) /DSQRT (XNN ( I ) ) 

CALL EIGENCB1.A1 , 6, 1 > 

DO 2 1=1,6 

2 Z1 (I )=DSQRT( A1 ( I ) ) 

DO 100 1=1, N 

DO 10 L= 1,6 

T2 (L)=BARN1 (-1 , 1 , 1 2787 , B2 ( L ) ) 

10 El (L )=BARN1 (-1 , 1, 12787, Z1 (L>) 

DO 89 LL= 1 , 6 
F=0 .ODO 
DO 85 'LLL=1 ,6 

85 F=F+(E1(LLL)*B1(LL,LLL>) 

89 T 1 (LL )=F 
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DO 15 L=l»6 
T 1 CL) =T 1 (L)+X CL) 

15 T2(L)=T2 (L)+XCL) 

CALL XNORMCT1 ,XM1 ) 

CALL XNORM (T2»XM2 ) 

DOT 1= ( XM ( 1 )*XM1 ( 1 > >+<XM(2)*XMl (2 ) ) + < XMC3)*XM1 (3) ) 
DOT 2= (XM ( 1 )*XM2( 1 ) ) + ( XM ( 2 ) * XM2 ( 2 ) ) + ( XM ( 3 ) *XM2 ( 3 ) ) 
FI 1 (I ) =D ARC OS ( DOT 1 ) 

F I 2 ( I ) =D AR CO S ( DO T 2 ) 

FI 1 (I )-FI 1 ( I >*57 .29578D0 
F I 2 { I )=FI2(I ) *57 • 29578D0 


100 

70 

CONTINUE 
WRITE (6,70 ) 

FORMAT ( 1H0, *V/E WILL 

NOW 

PRINT 

F 1 1 

RESULTS 

1 > 

71 

CALL FIND ( N » FI 1 ) 
WRITE (6*71 ) 

FORMAT (1 HO, 1 WE WILL 

NOW 

PR I NT 

F I 2 

RESULTS 

1 ) 


CALL FIND(N,FI2> 

RETURN 

END 

SUBROUTINE XN0RM(T1,M1) 

IMPLICIT REAL *8 (A-H,M»0-Z) 

DIMENSION T1(6)»M1(3> 

B 1 X =T 1 ( 1 ) 

B 1 Y=T 1(2) 

BIZ =T 1(3) 

B2X=T1 (4 ) 

B2Y=T1 <5 ) 

B2Z=T 1 C 6 ) 

XX=e2X-BlX 
Y=B2Y-B1Y 
Z— B2Z— B1 Z 

GX=(B1Y*B2Z )-(BlZ*B2Y ) 

QY= (B1Z*B2X )• — ( B1 X*B2Z ) 
GZ=(B1X*B2Y )-(BlY*B2X > 

MX = (GZ*Y )-(GY*Z) 

MY=(GX*Z )-(QZ*XX ) 

MZ= <GY*XX )~ (GX*Y > 
M=DSQRT(MX**2+MY**2+MZ**2) 

Ml ( 1 >=MX/M 
Ml (2>=MY/M 
Ml (3)=MZ/M 
RETURN 

END 

SUBROUTINE G AUSS ( /I X/ , S , AM * V, H ) 
IMPLICIT REAL*8 . ( A-H * O-Z ) 

K=H 

A=0 • OD 0 

DO 50 I = 1 * K 

CALL RANDUC IX, IY ,Y ) 

IX = I Y 
50 A=A+Y 

HO=H / 1 2 * 

H2=H/2* 

V=(S*( A-H2) ) /DSQRT (HO)+AM 
RETURN 

END 

SUBROUTINE RANDUC /IX/, I Y , YFL } 
IMPLICIT REAL*8 <A-H,0-Z> 

DATA JJJ5/1027/ 

IY=IX*JJJ5 
IF ( IY )5*6,6 

5 IY=I Y + 21 47483647+1 

6 YFL=I Y 

YFL=YFL*. 46566 2 3D-9 
RETURN 

END 

FUNCTION BARN 1(1*1 KEY * X FRN» SD ) 
IMPLICIT REAL*8 <A-H*0-Z) 


C 

C SD THE DESIRED STANDARD DEVIATION 

C AMEAN THE DESIRED MEAN 

C H THE POPULATION SIZE 


DATA AMEAN/O *D0 / 
DATA I HE RE /I 2787 / 
DATA H/36.D0/ 

I F ( I KEY ) 5 » 4 ♦ 4 

4 IHERE=I FRN 

5 IF ( I ) 6 *7 ,7 
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non non ro »- non nnnnn 


6 GALL GAUSS { IHERE»SD* AMEAN » VAL * H ) 

IFRN=IHERE 

GO TO 8 

7 CALL RANDU (IHERE, IFRN, VAL) 

IHERE=IFRN 

8 BARN 1 =VAL 
RETURN 

END 

SUBROUTINE EIGEN( AA, VALU.NR ,M) 

IMPLICIT REAL*8(A-H,0-Z) 

REALMS IND 

EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC MATRIX 

CRITICAL .NOTE EIGEN LIMITED . . .RANK MUST BE GE 2 OR LE 8 

NOT REALLY CAPABLE OF N X N MATRICES 


DIMENSION A<8,8),6<8,8 ) » VALU (8) ,DIAG (8 ) ,SUPERD(7 > , Q (7 ) , VALL (8 ) 
1,S(7),C(7),D(8),IND(8),U(8) , DUMMY! 94) , AA ( 64 ) 

EQUIVALENCE ( DI AG ( 1 ), DUMMY ( 1 )) , (SUPERDC 1 ) .DUMMY ( 9 ) ) * 

X (VALL C 1 ) »D( 1 ) , DUMMY { 16) ),<Q(1>,S(1) , DUMMY (24 ) > , 

X (B( 1 ) 1 ), DUMMY (31 ) ) , ( I ND ( I ) , U ( 1 > ) , ( 1 1 i MATCH) , 

X (TAU.8ETA), (P, PRODS 5, (T.SMALLD). ( ANORM . AN0RM2 ) 

CALCULATE NORM OF MATRIX 


N=NR 

ORMA = 0 .DO 
J =1 

DO 1 1=1 »N 

ORMA = ORMA+AA(J) 

J= J+N+ 1 
DO 2 1=1 .N 
NI=N*< 1-1 ) 

DO 2 J= 1 >N 
I J=NI+J 

A(J,I) = A A ( I J ) /ORMA 

3 AN0RM2=0.0D0 

4 DO 6 I =1 »N 

5 DO 6 J= 1 » N 

6 AN0RM2=AN0RM2+ACI , J)**2 

7 ANORM=DSQRT (AN0RM2) 

GENERATE IDENTITY MATRIX 

9 IF (M) 10', 45, 10 

10 DO 40 I =1 , N 
12 DO 40 J=1,N 
20 IF(I-J) 35, 25, 35 

25 8(1 ,J)=1.0D0 
30 GO' TO 40 
35 B(I , J)=O.ODO 
40 CONTINUE 

PERFORM ROTATIONS TO REDUCE MATRIX TO JACOBI FORM 

45 I EX I T = 1 
50 NN=N— 2 

52 IF (NN) 890, 170, 55 

55 DO 160 1=1, NN 
60 11=1+2 
65 DO 160 J=I I »N 
70 T1=A(I ,1+1 ) 

75 T 2= A ( I , J) 

80 GO TO 900 
90 DO 105 K=I , N 

95 T2=CUS*A(K,I+1 )+SUN£A(K, J) 

100. A(K , J)=CU$*A( K, J)-SUN*A(K,I + 1 ) 

105 A ( K , I + 1 )=T2 
1 10 DO 125 K=I , N 

115 T2=CUS*ACI+1 ,K)+SUN*A< J,K) 

120 At J,K)=CUS*A< J,K)-SUN*A(I+1,K) 

125 A(I+1»K) =T2 

128 IF (M) 130, 160, 130 

1'30 DO 150 K = 1 ,-N 

135 T2=CUS*B( K, 1+ 1 )+SUN*B(K, J ) 

140 B(K, J)=CUS*B(K, J)-SUN*B(K,I+1 ) 

150 B(K ,1+1 )=T2 
160 CONTINUE 
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non non non non on 


MOVE JACOBI FORM ELEMENTS ANO INITIALIZE EIGENVALUE BOUNDS 

170 DO 200 1=1 * N 
180 DI AG ( I ) = A ( I ,1 ) 

190 VALU ( I )=AN0RM 

200 VALL(I }=— ANORM 

210 DO 230 I = 2,N 

220 SUPERD (I— 1)=A(I— 1,1 ) 

230 GCI-1 )=<SUPERD(I-1 ) >**2 

DETERMINE SIGNS OF PRINCIPAL MINORS 

235 TAU = 0 ♦ OD 0 

240 1=1 

260 MATCH^O 

270 T2=0o0D0 

275 Tl=l .0D0 

277 DO 450 J=1,N 

280 P=DI AG( J )-TAU 

290 IF ( T2 ) 300 , 330, 300 

300 IF(T1) 310, 370, 310 

310 T=P*T1-Q ( J-l ) *T2 

320 GO TO 410 

330 IF(T1) 335, 350, 350 

335 T 1 =— 1 • OD 0 

340 T = -P 

345 GO TO 410 

350 T1=1.0DO 

355 T=P 

360 GO TO 410 

370 I F (O ( J— 1 ) > 380, 350, 380 
380 IF ( T2 ) 400, 390, 390 
390 T=— 1 m 0D0 
395 GO TO 410 
400 T = 1 # 0 DO 

COUNT AGREEMENTS IN SIGN 

410 IF(T1) 425, 420, 420 

420 I F ( T ) 440, 430, 430 
425 I F (T ) 430, 440, 440 
430 MATCH =M AT CH+1 
440 T 2=T 1 
450 T1=T 

ESTABLISH TIGHTER BOUNDS ON EIGENVALUES 
460 DO 530 K=1,N 

465 IF < K-MAT CH ) 470, 470, 520 

470 IF (TAU-VALL (K) ) 530, 530, 480 

480 VALL ( K ) = TA U 
490 GO TO 530 

520 IF (TAU~VALU(K) ) 525, 530, 530 

525 V ALU ( K ) =T AU 
530 CONTINUE 

540 I F { VALU ( I ) -VALL ( I ) -5 »0D~8 ) 570, 570, 550 

550 I F ( V ALU ( I ) ) 560, 580, 560 

560 I F ( DABS (VALL < 1 ) /VALU (I ) — 1 ♦ 0 DO ) —5* QD*~$ ) 570 , , 570 , 580 
570 1=1 + 1 

575 IF(I-N) 540, 540, 590 

580 T AU= ( VAL L C I ) +VALU (I ) ) /2 • 0D0 

585 GO TO 260 

JACOBI EIGENVECTORS BY ROTATIONAL T R I A NGULAR I ZAT I ON 

590 IF (M) 593, 890, 593 

593 I EX I T =2 

595 DO 610 1=1, N 

600 DO 610 J=1 , N 

610 A ( I , J ) =0 • 0D0 

615 DO 850 I =1 ,N 

620 IF ( I — 1 ) 625, 625, 621 

621 IF (VALU <1-1 >-VALU( I )-5,0D-7) 730, 730, 622 

622 IF (VALU ( I — 1 ) ) 623, 625, 623 

623 IF (DABS ( VALU ( I ) /VALU < I -1 ) -1 • 0D0 ) -5. OD-7 ) 730, 73 0r 9 6 25 

625 CUS=1 • 0D0 

628 SUM=0,0D0 

630 .DO 700 J=1 ,N 

635 IF(J-l) 680, 680, 640 
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640 GO TO 900 

650 S(J-1)=SUN 

660 C ( J-l )=CUS 

670 0( J-l )=T14CUS+T2$SUN 

680 T 1 = (DI AG ( J > -VALU { I ) ) *CUS-BET A*SUN 

690 T2=SUPERD<J) 

700 BETA=SUPERD < J )*CUS 

710 D ( N ) =T 1 

720 DO 725 J= 1 , N 

725 I NO ( J ) =0 » ODO 

730 SMALLD=ANORM 

735 DO 780 J=1,N 

740 IF ( IDINT( IND(J) )— 1 ) 750,780,780 

750 IF (DABS < SMALLD ) -DABS (D(J)))780. 760, 760 

760 SMALLD=D(J) 

770 NN=J 

780 CONTINUE 

790 I ND ( NN ) = 1 • ODO 

800 PRODS= 1 • 0 DO 

805 IF (NN-l) 810, 8S0 , 810 

810 DO 840 K=2,NN 

820 I I =NN+ 1— K 

830 A(II+1,I)=C(II)*PRDDS 

840 PR0DS=— PRODS'rS ( II) 

850 A ( 1 , I )=PRODS 
C 

C FORM MATRIX PRODUCT OF ROTATION MATRIX WITH JACOBI VECTOR MATRIX 

C 

855 DO 885 J=1,N 
860 DO 865 K=1 ,N 
865 U ( K )=A ( K» J) 

870 DO 8851 1=1, N 

875 A(I , J)=0.0D0 
880 DO 8852 K= 1 , N 

A(I ,J)=B< I ,K)*U(K)+A(I , J) 

8852 CONTINUE 
8851 CONTINUE 
885 CONTINUE 

DO’ 886 I — 1 » N 
NI=N* < I — 1 ) 

DO 886 J=1,N 
I J=NI+ J 

886 AA ( I J ) = A ( J , I ) 

890 CONTINUE 

DO 891 1=1 ,N 

891 VALU ( I ) = VALU(I)*ORMA 

RETURN 

C- 

C CALCULATE SINE AND COSINE OF ANGLE OF ROTATION 

C 

900 IF (T2 ) 910, 940, 910 
910 T=DSQRT ( T 1**2+T2**2 > 

920 CUS=T1/T 
925 SUN=T2/T 

930 GO TO (90,650), IEXIT 
940 GO TO (160,910), IEXIT 
RETURN 
END 
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APPENDIX D 


SAMPLES OF OUTPUTS OF PROGRAMS 
D.l Shock Progam Output 

The output, as it appears in printed .form (the upper portion on the next two 
pages; the lower portion appears on the two pages after that). The sample shown 
is for a simulated shock example and is fully described in Section IV-B of the 
report. A real shock sample output would be almost identical except that also 
printed out in the upper portion would be the Y array, and XMEAN does not usually 
replace XSTART as in the simulated cases, but it can. 

D.2 Cone Program Output 

A sample of this output (appears on the single page following the Shock 
Program output sample). Section IV-B also fully describes this printed output 
for simulated shocks. It is identical in appearance for real shock cases. 
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GPAM— QUTPtJT 


SIMULATED SHOCK BEING PROCESSED 


THE NUMBER OF DIFFERENT CASES TO BE PROCESSED IS 


T HE TOTAL NUMBER Op DATA POINTS, N, IS 140 



mi:- in i — Diu — w a -j 

THE INPUT NN1 WAS 



20 

— U . '_IU u 

20 

B.JJU 

20 

. O J V 

10 

INPUT XST ART WAS 

4.000 

5.000 

~ l . 000 

3.500 

the corresponding g 

VALUES ARE 

T J V* 




N I N 2 

7.260 13.861 

wx 

75. 633 

N X 

0. 953 

N f 

0.222 

N £ 

0.207 


THE NUMBER Op ITER AT I ONS *L . I S 



THE BEST 

ESTIMATE INDEPENDENT 

»— > 1 -r 

PARAMETER 

MATRIX is 


B 1 X ' 
3. 99 8 

OAT 

5.095 

D 14 

-0.977 

B2X 
3. 660 

02 T 
9.475 

02 c. 

-3.829 


3. 99 8 

5 .095 

— Hi DTQ 

-0. S77 
— 

3. 8 6 0 
3 i 

9.475 

Q i6no 

—3 • 8 29 
A T 


4,02 1 
- 4,026 

4.026 


4.026 


4. 026 


4.026 


4. 026 


4.026 


5 * 0 57 
—5 »-0-S5- 

5 . 056 


5.056 


5 o 0 56 


5.056 


5.056 


5.056 


- 1 .033 
Hr^-0-25- 
~ 1 ♦ 024 


-1.024 


-1.024 


- t. 024 


- 1.024 


- I . 024 


3. 754 
3. 723 


3 © 723 


3. 723 


3.723 


3. 723 


3. 723 


9.649 

- 9 -y 6 z* 0 - 

9.615 



9.614 


-3.069 
-3- » 1 9 4- 
-3.2 13 


5 


215 


5 


-3.215 
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82 Y 

B2Z 

WY 

WZ Ni 


N2 

WX 

L « 1 0 0 
10 

1 . 300 

10 

"1 0*000 

10 

— 1'OTrOOt) O. rOU 

10 10 

1 * 

0 00 1 O" 

10 

• Uuu 

10 

9.000 

-3.000 

1 0.000 

20.000 


* 



— — 

13*361 








- 

WY — 

4 * a^5 

w z ~ 

17.983 

N 1 

6.890 

w a 

13.4-86 71.404 MEAjg VALUES 


~ o /tj r • 

4*255 

17*983 


8707.94392=LOSS 0 

Z= 

0*2197 

^ A_.. 


0 . 930 

4 .070 

g 2~»~0~0~3 — 

20.107 


— 2 to #”7 1 0(iJ“LUSo T“ 
163.04689=LDSS 2 

< pCO-i-n.Cf' "3 

2.— 

z= 

v 4 1 JO 4 

0. 0931 
r\ r>F> n o 


3“* 515 
3*499 

21*296 


13 r * OOODL-Luoo O 

1 36 • 976 95 —LOSS 4 

' z 

z= 

0.048 7 


• ■ “UTUvS* 

3*504 

21 *323 
21.324 


X *5 O • V # v f w L LJ v v v* 

1 36 *97 673=L0SS 6 

Z— 

7 — 

v • 04-3 y 

0.041 3 

/> a a n r> — 


3‘»~50-4- 

3 *50 4 

4, 1 . o Z 4 
21 .324 


iJutynj/ o^Lub'ts r - 

136 .97673— LOSS 8 

z= 

V . 040 £ 

0 . 0397 


3'*50 4^- 

3*504 

— 2 1 « 324 — 
21.324 


1 36 • 9 f o r“3-tu^u 7 

136 .9767 3= LOSS 10 

.'t r- __r\rr .ei \ 

Z — 

z= 

U ♦ U - W 4 

0.0393 


3<5v4' 

3*504 

2 1 *Oc4 

21 .324 


i ou * v 7 a r^-i_ao o ri z_ — 

136 .97 67 3= LOSS 1 2 Z= 

U • U Jy c 
0.0392 


3 * 50 4 

3 *504 
3-*50-4- 

21. 324 — 

21.324 
— 2’l~r3-24-- 


0'-l_U0C5 1 -3 

136 *9767 3= LOSS 14 
~i^-^7673^t-0S5 -T5- 

— 

U . U^>9 tL 

0.0392 
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Tf-Hr — — EST- i HATE — DEP EMDGN T-- PARAMETER MAtTRJ-X-- 1-S 

N 1 N 2 hX NX NY ' NZ 

MgAN»S G f S 3 6, 018 3 6 , 737 0.003 0.3 6 4 0 . 3 4 7 — 


18. 943 
7.332 -- 
7*510 
7* 0 00 


38.918 
14.048 
14. 647 


36. 
1 07. 
77, 

■ - 7-I-. 
71. 
-7*1. 
71 . 
71'* 
71 . 

— -jH-t* 
71 . 
71. 
71. 
71* 
71* 
71. 


737 

703 

973 

f>tr ;■> 

097 
-0-72 ■ 
071 
071 
C71 
>0 74 - 


0.903 

0.977 

0.967 

0*9 5 9 


0.264 
0. 1 15 
0. 148 
0* 170 


0.347 
0 * 130 
0.208 
0*225 


6.931 
6.932- 

6.932 

- 6,932- 
6*932 
6 *9 32- 
6.932 

-*• 6r. 9 32 
6.932 

- 6,932- 
6*932 
6.932 


13.457 

-13,459— 

13.438 

13*468 

13.458 


13.458 
13, 45a * - 
13.458 
V3* 458 
13. 456 
13.458 


0.958 
— 0,958* 
' 0.958 
- 0*458 
0.958 
0*056 


0. 174 
4>. 174- 
,0. 174 
0.1-74 - 
0.174 

" & y l - 7 ‘4 

0. 174 

-0-.-1-74 

O. 1 74 
0,-174- 
0*174 
0.174 


0 . 229 
- O • 229 
0 .229 

- 0.229 
0.229 



0.229 

O * 229- 

O .229 

- 0 . 229 
0.229 
0 . 229 


C71 
C-7-1 - 
C7 1 
07 Ir- 
CH 
C71 


0.953 

*-0.-950- 


0. 958 
0-058- 
0.958 
0.953 


-'3-§~-T+IE-COVAR-f ANCE -MATRIX OF F f* NA i_ ESTj MATE »' 'FS' 

0.005B5 0.00016 0.00018 0.00173 0.00012 

O-rOO-O-t-B -O-.-0-L-Q68 0-. 00021 0'S 00-07-1 OtOO&SA- 

0.00016 -0.00021 0 .00556 -0.00057 0.00104 

0 .0 O T 7 3 “- 0 - . - 000 7 1 — O.OOC57 0.02461 »0. - 002 - 9? - 

0.00012 0.00584 0.00104 -0.00207 0.07410 

O - .OO - lO - l - 6~r€H>284 0.-00 7-9-7 <3>-O0 2^1 -0^-01-690 - 

0.02443 0.03608 0.00433 -0.18015 -0.29256 

—0 ~ r 02-72 4 O-s-0 1^494 0 tt02-770 0 *15 06 0 0t0O4:>2 - 


-0.00181 0.02443 

— ©-S-002Q4 0.0 3608 “ 

0.00797 0.00433 

0.00291 ^0. 1 8015 

-0.01690 -0.29265 

— O-y-05650 0-yl-50-t3- - 

-0.15013 5.74353 

—0 - 3 3-18^ 2^0-1-930-- 


T ' l !C "GO NT 1 R ACT CD TOflH 3F DEfM ' VA - T - i-V E - MA T3 I X > A > tO 


1,000 

*■ * O.-0 

0. 0 

0.0 

0. 0 

o, c 

0.0 

0.0 

0 . 0 

i'i“6O0- 

0.0 

i-N a : . 

— O'. O' ~ 

t . 000 

. ' A.. A 

0 « O' 
0. 0 

o.o 

0.0 

0.0 

0.0 

0“5"0 

0.0 

0"? - C - 

0.0 


It "000 — 

0*0 

0*0 

0,0 

0*0 

ftrfl 

0. 0 

— a- 

0*0 

8r-e 

0, 0 

a-r-e 

1 * 000 

Bt-0 

0*0 



0 *0 

~~ 


0.0 0.0 0 . 0 0.0 0.0 0.0 1.000 

— Os-6 O-sO OsO ChrO OsO- OsO Os"0 — - 

-10.618 2.474 12.937 9.962 0.145 -5.866 0.445 


20,-61-4 

— 2-*-2-26 

— 25-*-009— 

IS* 340 

1-^608 

l"t-,"6~06 

0-4^39- 

33.718 

47*999 - 

104, 491 

- 17,369 

24. 724 

53*824 

1 *424 

AAV E= 8.7 34 

A8E= 

3.-0 88 

AVE ,ABE= 

8.772 
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DP 

7t 62 4 


0»12 6 fr e- t5HM!_I TY M 


7,624 

27,292 

13,634 

10,007 

10*003 

10.003 

10.003 

10.003 
-M>irOO-3- 

10.003 
10.003 
10.003 
10.003 
10.003 
10.003 


0 . 12680=0 UA LI TY 0 

0 . 70 8^=aWL-m ir 

0 • 92663 = QUALI TY 2 

1 * 009 5 3 — QUA LI TY =h 

1 .01097=QUALI TY' 4 

-1-^01096 —QUA t_J-T ¥ 5" 

1.01 098=G UA LI TY 6 

' - 1,01 098=0 UAbf-TY 7- 

1,01098=QUALITY 8 

1 .01 09 8- ft - UAiii TV 9 

1 *01098=GU\LI TY 10 

1 .0 1098=0 UA LI TY 12 

1^ 01-09 B = 3 UA bl-T-Y- -1-3- 

1 .01098=QUALI TY 14 
1 »0109B=QUALI TY 15 


-0 ,0 2724 

0.0 1494 

0 .02770 

O. TO Trt rO 

0.00422 

0.-33189 - - 

2.01930 

7 . 10262- 


0*0 




0.0 

0.0 




V « u 
0,0 

U • 

0,0 

“ i'.-OOO- * 

-0.723 


• 





TT4rT “ " 

3.099 
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THE— INPUTS 

THE INPUT 

X VALUES-* 
S1G VALUES 

frfrE 

WERE 0.35 

5.00 1- 

0.50 0 

5.50 

.35 0.60 

9.00 3.00 

1.10 1.30 

TTTE — INPUT 

THE INPUT 

VALUE — FOR 
NN VALUES 

IN WAS 3000 
WERE 20 

20 

20 10 

10 10 

THE INPUT 

VALUES FOR 

B 1 ' WERE 




Q • 0 C 5 8 5 
0.00016 

0.00016“ 
0.01 063 

0.0001 8 

-0.00021 - 

0.C01 73 
0,00071 

0.00012 

0.00884 

— 0 #00 1 81 
0. '00 284 

A ..7\.A.7.A7 

0.00018 

0.00 173 

— 0.000 2T 
-0.00071 

0 • 00 556 — 

-0.00057 

0*00057 
C. 02461 

• 0.001 0"4 
-0.00297 

O • v O fy/ 

0.00291 

0 • 0 C C 1 2 
-0.00181 

0 *uOd8A 
0.00284 

0.001 04 — 

0.00797 

0*3 029 7 
0.00291 

0 . C"74T 'V' ‘ 
-0. C1690 

■■ "0*0 1 uyu 
0.05350 


WE WILL NOW PRINT FI1 RESULTS 


ISIG1= 150 1 S I G2= 30 ISIG3= IS 


THE lb VALUE WAS 0.700376710 01 


THE 30 VALUE WAS 0.6501Q9930 01 


THE 150 VALUE WAS 0.51747948D 01 


WE WILL NOW PRINT FI2 RESULTS 


I S I G 1 = 150 ISIG2= 30 ISlG3= IS 


THE 15 VALUE WAS 0.14611719D 02 


THE 30 .VALUE WAS 0.13564957D 02 


THE 150 VALUE WAS 0.106346170 02 
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